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COMPUTER PROGRAM FOR CALCULATING VELOCITIES AND STREAMLINES 
ON A BLADE-TO-BLADE STREAM SURFACE OF A TURBOMACHINE 

by Theodore Katsanis 
Lewis Research Center 

SUMMARY 

A FORTRAN IV computer program was written that gives the solution of the two- 
dimensional, subsonic, compressible (or incompressible), nonviscous flow problem for a 
rotating or stationary circular cascade of blades on a blade -to -blade surface of revolu- 
tion. The flow may be axial, radial, or mixed. There may be a change in stream chan- 
nel thickness in the through-flow direction. 

The computer program requires the basic cascade geometry, the meridional stream 
channel coordinates, fluid total conditions, weight flow, and inlet and outlet flow angles. 
The output includes streamline coordinates, velocity magnitude and direction throughout 
the passage, and the blade surface velocities. 

The method is based on the stream function with the solution of the simultaneous, 
nonlinear, finite-difference equations being obtained by two major levels of iteration. 

The inner iteration consists of the solution of simultaneous linear equations by successive 
overrelaxation, using an estimated optimum overrelaxation factor. The outer iteration 
then changes the coefficients of the simultaneous equations to compensate for compres- 
sibility. 

This report includes the FORTRAN IV computer program with an explanation of the 
equations involved, the method of solution, and the calculation of the velocities. Numeri- 
cal examples have been included to illustrate the use of the program, and to show the re- 
sults which are obtained. 


INTRODUCTION 

In the design of blade rows for turbines or compressors, it is desirable to obtain the 
velocity distribution through the passage and particularly over the blade surfaces. The 
trend to highly loaded blading results in more widely spaced blades with less of the pas- 



sage being within a guided channel between the blades. The velocity distribution is read- 
ily obtained within the guided channel by stream filament techniques. 

For the unguided portion of the passage, finite-difference methods have been used. 
Stanitz (refs. 1 and 2) obtained finite-difference solutions for compressible flow through 
turbomachines, without the use of a computer. Kramer (refs. 3 and 4) has obtained 
finite-difference solutions for incompressible flow through centrifugal pumps, using a 
computer for the solution of the finite -difference equations for the stream function. More 
recently a program has been written to perform in addition to the solution of the finite- 
difference equations, the calculation of the coefficients, and the differentiation of the 
stream function to obtain the velocities for incompressible flow through an axial blade 
row (ref. 5). 

To extend this technique, a computer program has been written to obtain a numerical 
solution for ideal, subsonic, compressible (or incompressible) flow for either an axial, 
radial, or mixed flow turbomachine blade row which may be fixed or rotating. The stream 
function used here is a function of meridional streamline distance and angular coordinate, 
whereas the previously mentioned references all used either radius or axial coordinates 
instead of the meridional streamline distance. Also, the finite -difference equation has 
been given in a simpler form. The input required consists of the basic geometry coor- 
dinates, fluid total conditions, weight flow, and inlet and outlet flow angles. The output 
includes velocity magnitude and direction through the passage, blade surface velocities, 
and streamline coordinates. 

This report includes the FORTRAN IV computer program that was developed, with an 
explanation of the equations involved and the method of solution. A radial gas turbine 
rotor and an axial turbine stator have been analyzed to illustrate the use of the program, 
and these results are compared with results obtained by other methods. 

This report is organized so that the engineer desiring to use this program needs to 
read only the sections MATH ANALYSIS, NUMERICAL EXAMPLE, and DESCRIPTION OF 
INPUT AND OUTPUT. The necessary information of interest to a programmer is con- 
tained in the sections DESCRIPTION OF INPUT AND OUTPUT and PROGRAM PROCE- 
DURE. 


SYMBOLS 


A 

a 0> a i> a 2’ a 3 
a 4’ a 12’ a 34 


coefficient matrix, eq. (A7) 
coefficients in eq. (A 2) 


a ij 


typical element of matrix A 


2 



normal stream channel thickness, m 
quantities in eq. (A 2) 

specific heat at constant pressure, J/(kg)(°K) 
spacing between adjacent points, eq. (Al), see fig. 18 

, eq. (A 7) 

coefficient matrix of eq. (A8) when 0 = 1 
meridional streamline distance, see fig. 2 
number of unknown mesh points 
gas constant, J/(kg)(°K) 
radius from axis of rotation, m 
angular blade spacing, rad 
temperature, °K 
stream function 

discrete approximation to stream function at n mesh points, 

m^ iterate of u, 

absolute fluid velocity, m/sec 
fluid velocity relative to blade, m/sec 

mass flow per blade flowing through stream sheet, kg/sec 
axial coordinate, m 

angle between meridional streamline and axis, rad, see fig. 

angle between relative velocity vector and meridional plane, 
see fig. 1 




y specific heat ratio 

7] outer normal to region 

9 relative angular coordinate, rad, see fig. 1 

2 

A. prerotation (rV 0 ) , m /sec 

y in 

O 

p density, kg/m 

p( ) spectral radius of matrix 

f2 overrelaxation factor, eq. (A8) 

a) rotational speed, rad/sec 

Subscripts : 

cr critical velocity 

giv given 

i dummy variable 

in inlet or upstream 

j dummy variable 

1 lower surface of blade 

m component in direction of meridional streamline 

out outlet or downstream 

r radial component 

u upper surface of blade 

z axial component 

9 tangential component 

0,1, 2, 3, 4 quantities at these locations in finite difference expression, fig. 18 
Superscripts: 

T transpose of vector or matrix 

’ absolute stagnation condition 

,T relative stagnation condition 


4 


MATHEMATICAL ANALYSIS 


r 


It is desired to determine the flow distribution through a stationary or rotating cas- 
cade of blades on a blade- to- blade surface. The following simplifying assumptions are 
used in deriving the equations and in obtaining a solution: 

(1) The flow is steady relative to the blade. 

(2) The fluid is a perfect gas or is incompressible. 

(3) The fluid is nonviscous. 

(4) There is no loss of energy. 

(5) The flow is absolutely ir rotational. 

(6) The blade-to-blade surface is a surface of revolution. 

(7) The velocity component normal to the blade-to-blade surface is zero. 

(8) The stagnation temperature is uniform across the inlet. 

(9) The velocity magnitude and direction is uniform across the upstream and across 
the downstream boundaries. 

(10) The relative velocity is subsonic everywhere. 

The flow may be axial, radial, or mixed and there may be a variation in the stream chan- 
nel thickness b in the through-flow direction. 

The coordinate system is shown in figure 1. Since the variables r and z are not 
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Blade-to- 
blade sur- 
face 

Figure 2. - Blade-to-blade surface of revolution. 

independent on the stream surface, one variable can be eliminated. It is better, how- 
ever, to use the meridional streamline distance m as an independent variable (see 
fig. 2). Then, m and 0 are the two basic independent variables. A stream channel is 
defined by specifying a stream channel thickness b as shown in figure 3. 

For the mathematical formulation of the problem the stream function is used. The 
stream function u satisfies 

_L 13 \ d 2u - Jl i l£ 1^ + sin a _ J_ dfop) _lH_ = 2b P^ sin a (!) 

2 ,2 2 p 30 38 r bp 9m 3m w 

r 00 0m r ^ u ' - 1 

This may be obtained from equation 12(9) of reference 6 by letting u = -i^/w, where \(/ 
is the stream function as defined in reference 6. The stream function u has the value 0 
on one blade and 1 on the other. Also, the derivatives of the stream function satisfy 

iH_=- b P w„ ( 2 ) 

0m w i 

1H = W (3) I 

00 w m | 
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1 1 

Figure 3. - Flow in a mixed flow stream channel. 



Figure 4. - Finite flow region for a radial turbine. 


For the solution of equation (1), a finite region is considered (as indicated in fig. 4) 
with the condition that the flow along AB is the same as along HG, and the flow along CD 
is the same as along FE. Also, it is assumed that AH is sufficiently far upstream so that 
the flow is uniform along this boundary, and that the flow angle /3- n is known. Similarly, 
it is assumed that the flow is uniform along DE, and that the flow angle /3 QU £ is known. 
For an actual blade row, /3 ^ may usually be determined by means of experimentally de- 
termined rules. Specifying /3 ^ along DE is mathematically equivalent to specifying 
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the location of the stagnation point on the trailing edge of the blade. 

Since equation (1) is elliptic for subsonic flow, boundary conditions for the entire 
boundary ABCDEFGH are required. Along BC, u = 0; along FG. u = 1. Along AB, GH, 
CD, and EF, a periodic condition exists; that is, the value of u along HG and FE is ex- 
actly 1 greater than it is along AB and CD. Along AH and DE, du/dri is known, where 77 
is in the direction of the outer normal. From equations (2) and (3), since W^/W m =tan /3, 

tan/3 (4) 

3m rdO 


Along AH and DE, 


so that 


3u _ u(H) - u(A) _ _1 
30 s s 



tan /3. 


in 


sr 


in 


along AH 



ton ^out 
sr out 


along DE 


(5) 

( 6 ) 


These are the boundary conditions required to determine a solution to equation (1). The 
method used for the numerical solution of equation (1) is described in appendix A. 

After computing a numerical solution to equation (1) in a given flow region, the veloc- 
ity at any point can be computed from equations (2) and (3) by using numerical differen- 
tiation. The streamlines are located by the contours of equal stream -function values. 


NUMERICAL EXAMPLES 

To illustrate the use of the program and to show the type of results which can be ob- 
tained, two numerical examples are given. The first example is a radial inflow turbine, 
and the other is an axial turbine stator. 
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Radial Inflow Turbine Rotor 


.f 


The turbine profile is shown in figure 5 with the mean streamlines and mean stream- 
sheet thickness as calculated by the quasi -orthogonal method (ref. 7). The blade-to- 
blade shape in 9 and m coordinates is shown in figure 4. This particular rotor had 
splitter blades as indicated in figure 5. There were 11 complete blades and 11 splitter 
blades. The program cannot handle the case where the blades are not all identical. 
However, two solutions can be obtained, one based on 11 blades and one based on 
22 blades. The solution with 11 blades should be reasonable for the region beyond the 
splitter blades, and the solution with 22 blades should be reasonable for the region with 
the splitter blades. Note that this technique would not work for a compressor with split- 
ter blades, since the percentage of flow on each side of the splitter blade would not be 
known. 



Figure 5. - Flub-shroud profile with streamlines used for blade-to-blade analysis of a 
radial turbine. 


The input for the case with 11 blades is given in table I. Because of a high local 
velocity near the trailing-edge radius, the program did not converge at every point (see 
error condition (6) (p. 36) for further discussion). However, the solution failed to con- 
verge at only one point, so that the rest of the velocities should be valid. 

The results are plotted in figure 6. There is also shown in this figure a solution ob- 
tained by the quasi -orthogonal method of reference 8. To make the results comparable, 
the quasi -orthogonal solution was obtained for zero loss. There is close agreement on 
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TABLE I. - ELEVEN-BLADE RADIAL INFLOW TURBINE ROTOR AND COMPUTER INPUTS 


GAM 

1 -66670C0 
CHORD 

G.684400QE-01 

PI 

0.648CCC0F-G3 
MXB I MX EC MX KR 

5 29 50 

MU ARRAY 

-0 

C.53940Q0E-01 
XSPU APR 

-0 

-0.962CCC0E-0I 
ML ARRAY 

-C 

C #52 940006-01 
XSPL APR 

-C 

-0 -1452COO 

HP ARRAY 

-C.7620000F-0? 

0 .46 23OCOE-0 1 
RMS P ARR 

C.E4C1000E-01 

O.44350C0E-01 
0ESP APR 

C.S70QOCOE-0? 

0.132OC00E-07 
BLDATA KULAK I Ei 
1 C 


AP 

2C 8.2CCCC 
STGR 

-C.539CCCC 

ALLI 

2.CCCCCCC 
3 I NLSF KLSP KRSP 
15 11 I! 16 

TIP 

1OE3.CC00 

BETAI 

54.2CCCCC 

ALLI 

■2.CCC0CCC 
KRL MKT 
11 5 

RFOI P 
0 .3956600 
eETAC 

-61.50CCC0 

PC 

0.753CCC0E-03 

V.TFL 

0.12580G0E-C2 

ALUO 

-63.80CCC0 

OMEGA 

4C3C.0C0C 

ALLO 

-63.8CCCCC 

Vi 

-C 


C.86CCCCCE-02 
0.6197CCCE-C1 - 

C. 16 CCCC0E-C1 
C 

C •2350000F-C I 

C.2903CC0E-C1 

C.3428CCCE-C1 

C.3S54CCCE-C1 

0. 4623000 E-0 1 

8 Y 

0.1220C0CE-CI 

-0.2791CCC 

0 • 1560CCOE-C1 
0.5427COO 

C.109OOCOE-C1 

C.2C9O0C0E-01 

0 .2 1 50CCCE- C 1 

C. 1600C00E-C 1 

-0. 11600C06— 01 

C .86C0CCCE-C2 
0.6157CCCE-C1 - 

C.16CC0C0E-01 

0 

0.2350CCCE-C1 

0.2903CC0E-01 

0.3428CCCE— 01 

C.3554CC0E-C1 

0 • 46 2 30 00 E-C 1 

»Y 

-0.12200C0E-C1 

-0.33470CC 

0.1560GCOE-01 

0.5842CC0 

— C . 19CCCC0E— C 1 

— 0.2C7CCCCE— 01 

-0.2220CC0E-C1 

-0.2840CC0E-CI 

-0 • 5580CC06— C 1 

0 

0.5394CCGE-CI 

0.86C00006-0? 

0.6197CC0E-CI 

0 . 1600000 E-C 1 
0.6 8440 00F-01 

C -2350CCCE— 01 
0.735 4000E-01 

0.2SC3C0CE-C1 

0.8116CCCE-CI 

C.3428C0CE-CI 

0.8878CC0E-CI 

0.39540 COE— 01 
0. 970 00 00 E-0 1 

! Y 

G.7645CCCF-CI 

0.4295CCCE-C1 

0.68CC00OE-01 
0.413 1C00E— 01 

0 .6 1C 3000 E-C 1 
0 .4C05CCOE— 01 

C.5471CCCE-01 
0.396 4000E-01 

C «5C£ SCCCE-C 1 
0 . 394 4CGCE-C 1 

C.48CECCCE-C1 
C.3940CCOE— 01 

0 -4602000 E-0 1 
0.3940000E-01 

! Y 

C.S6CC0CCE-C3 
0 • 14 300C0E-02 
5PRT STPFN SLCPD 

0.1C3CCC0E-02 
C .153 CO COE- 02 
i APPRT INTVEL 

C.1C90000F-C2 
0. 1620000 E-02 
SIP Y El 

C.l 14CCCOE-02 
0. 1670000E-02 

0.1 16CCC0E-C2 
C.16SCC0CE-02 

C . 1 16CCC0E-C2 
0.1 7 COG CQE-C2 

0. 124000CE-02 
0. 17000 COE— 02 



0 2 4 6 8 


Meridional streamline distance from leading edge, cm 
Figure 6. - Blade surface velocities for a radial turbine with splitter blade. 
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Figure 7. - Axial stator blade for numerical example. 


TABLE H. - AXIAL STATOR COMPUTER INPUT 


GAM 

AR 

TI P 

RH01 P 

*T FL 

OMEGA 


1 .AOUOOOC 

287.0530C 

2 88. 15000 

1.2250000 

0 .3146000 

0 


1.8500000 

CHORD 

STGR 

BETAI 

BET AC 





0.4265000E-0 1 - 

0. 1116150 

0 

-6 7 • OuOOoO 





R i 

ALU I 

ALL! 

RO 

ALUO 


ALLU 


0 .38 100 00 E- 0 2 

28.3C0C00 

-14.20C0QU 

U.889000UE-03 

-72 .400000 

-56 

. 1C0CC0 


MX B I MX 00 MX NBB 1 

NUSP NLSP NR SP 

NBL NUT 






15 32 47 20 

7 6 2 

50 5 






MU ARRAY 








-0 

0.8575000E-02 

0.1 71 5 QCOE-Ol 

0 .257 2500 E-0 1 

0 .3430000E-01 

0. 

38588CCE-01 

-0 

XSPU ARRAY 








-0 

0.1 7£9GCGE-C1 

O.15380C0E-G1 

-0.5310000E-02 

—0 .4654000 E-G 1 

-c . 

74OCCC0E-O1 

-0 

ML ARRAY 








-0 

C.85750CCE- C2 

0. 1 71 50 COE— G1 

0.2572500E-01 

0. 34300 00E-01 

-0 



X SPL ARRAY 








-0 

0.15620 COE- 01 

-0.2 854C00E-01 

-0.5070000 £-0 1 

-0 .B250000E-01 

-0 




MR ARK AY 

-0.500GC0CE-C1 C.IGGGGOO 


RMSP ARRAY 

Q.3302G0C C.23C200C 

bESP ARRAY 

0.101600C C.10160CC 

BLOATA NULAKI ERPRT STRFN SLCRO ARPRT I NTVEL SURVEL 
1C312013 
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TABLE m. - AXIAL STATOR COMPUTED SURFACE VELOCITIES 
(a) Surface velocities based on axial components 



* 


UPPER SURFACE 


* 


LOWER SURFACE 


* 

H 

* 

VELOCITY 

ANGLE (DEG) 

SLRF. LENGTH 

RHC*W 

* 

V6L0C 

TY ANGLE (DEG ) 

SURF. LENGTH 

RHC*W 

* 

0 

* 

C 

50.00 

0 

-6.3313 

* 

0 

-50. CO 

0 

-8.9717 

* 

0 .2509E-C2 

* 

57.625 

27.31 

C.4405E-02 

114.73 

* 

76.541 

-2C.67 

0.43726-02 

91.409 

* 

0 .50186-02 

* 

111.01 

21.66 

C. 7165E— 02 

128 .86 

* 

64.625 

-13.60 

C.6959E-C2 

77.745 

* 

0 .7526E— C 2 

* 

128. 9C 

14.63 

C. 98076-02 

146.82 

♦ 

60.552 

- 16.21 

0.9547E-C2 

73.054 

* 

0 .10046-01 

* 

148.82 

6.2C 

0. 1236 £-0 1 

165.37 

♦ 

59.431 

-21.52 

0.12196-01 

71.698 

* 

0 .12546-0 1 

♦ 

168.62 

-3.35 

C. 14876-01 

182.13 

* 

60.567 

-26.75 

0.1452E-01 

73.491 

* 

0. 15056-C1 

* 

185.87 

-13.61 

0. 1741 fc-01 

195.23 

* 

64.056 

-31.55 

0. 1 777E-01 

77.085 

* 

0 . 17 56E-C 1 

* 

195.22 

-23. 76 

0.20066-01 

204.35 

* 

68.661 

-35.51 

0.2075E-01 

82.4 08 

* 

0 .20076-01 

* 

2C6.91 

-33.36 

0.2292 6-01 

209 .18 

* 

74.60C 

-25.67 

0.23876-01 

89.205 

* 

0 .22506- C 1 

* 

211.94 

-41.92 

0.26096-01 

212.17 

* 

82.C35 

-43.47 

0.271 6E— 01 

97.602 

* 

0.25096-01 

* 

211.49 

-49.21 

G.2969E-01 

211.90 

* 

91.115 

-46.74 

0. 3063E-01 

107.66 

* 

0 .27606-01 

* 

205.15 

-54.84 

0.3379E-01 

208.10 

* 

1C1.5C 

-45.74 

0.342 86-01 

119.31 

* 

0.30116-01 

* 

156.59 

-58.36 

0. 3837 fc-01 

203.94 

* 

114. 6C 

-52.57 

0.381 4E -01 

132.55 

* 

0.32616-01 

* 

195.85 

-60. 50 

0. 43316-01 

202.14 

* 

129.22 

-55.22 

0. 42236-01 

147.24 

* 

0.35126— Cl 

* 

202.02 

-61.95 

0.48526-01 

206.15 

* 

146 .40 

-57.62 

0. 46566- Cl 

163.20 

* 

0 .3763E-01 

* 

205.68 

-66. 76 

0. 5428E-C 1 

210 .84 

♦ 

164.56 

-56.75 

0.51086-01 

179.18 

* 

C. 401 46- C 1 

* 

192.93 

-71.68 

0.61506-01 

200 .16 

* 

170. C2 

-56.49 

0. 55666—01 

183.24 

* 

0 .42656-01 

* 

0 

-50. OC 

C.6989E— 01 

51.236 

♦ 

0 

50. CO 

0.583 OE -01 

-47.293 

* 




(b) Surface velocities based on 

tangential components 






UPPER SURFACE 




LOWER SURFACE 


M 

VELOC 

1 TY ANGLE (DEG) 

RHC*W 

M 

VELGLJ 

[TY 

ANGLE (DEG ) 

RHC + W 

-o 

10.858 

50.00 

13.2 94 

0 .615QE-C3 

51.65C 


-56. 57 

o2 • 545 

0 .6157E-C2 

78.376 

56. 97 

93.484 

0 .4764E-C 2 

38.457 


-13.25 

46.810 

0 .25866- C 2 

55.856 

25.04 

112.87 

0.1139 fc-C 1 

60.251 


-23.55 

72.656 

0 .19066-C 1 

204.51 

-29. 55 

2C7.71 

0.1539 E-C 1 

64.57C 


-31. C7 

78.146 

0 .22016-0 1 

2C5.76 

-40. 05 

210.89 

0.1B5OE-C 1 

7C.9C3 


-36. 17 

84.983 

0 .2419E-C 1 

212. 11 

-46. 76 

212.26 

0.2114 6- Cl 

77.737 


-40. U5 

92.762 

C. 25986-01 

205.67 

-51. 5C 

21C.95 

0.234 76- Cl 

65.262 


-43.17 

101.20 

0.2753 6- Cl 

205.28 

-54. 73 

208.18 

0.2558 E-C ) 

53.344 


-45. 76 

110.09 

0 .28S4E-0 1 

2GC.81 

-56. 53 

205.37 

0.27 536- 0 1 

101.76 


-4 7. 9 5 

119.16 

0.3025E-C 1 

197.72 

-58.52 

203.37 

0. 29336-C1 

11C. 64 


-45. 97 

128.49 

0.3149E-C1 

196.38 

-59.68 

2C2.49 

0 .31026- 0 J 

119. 5C 


-51. 75 

137.93 

0.3268E-C1 

155.76 

-60.54 

202.08 

0 . 326 16- C ) 

129.46 


-53.35 

147.36 

0.2384E-C1 

197. C2 

-61.16 

202.91 

0.341 1E-C 1 

135.13 


-54. 81 

156.55 

0 .345 76- C 1 

2C1. 14 

-61.82 

2C5.59 

0 .35546— C 1 

145.55 


-55. 56 

166.05 

0 .3605E-C 1 

206.82 

-63.34 

2C9.12 

0 .36526- C 1 

155. 6C 


-56. 6 C 

1 74. 88 

0.37046-C 1 

2C8.57 

-65.38 

21C.42 

0. 28286-C1 

165.65 


-56. 83 

1 82 . 95 

0. 37556-C1 

206.22 

-67.52 

2 1C. 04 

0.3964E— C I 

178.75 


-56. 67 

190.00 

0 . 38766-C1 

204.55 

-65. 51 

207. 76 

0.4102E-C 1 

195.67 


-56. 1C 

2C4.64 

0. 39 5 IE- C 1 

199.01 

-70.85 

204.21 

0 .42506- C 1 

60.C78 


56.49 

95.402 

0 .4021E-C 1 

152.51 

-71. 74 

199.88 






0.4088E-C 1 

167.24 

-72.25 

196.21 






0 .4154E-C 1 

162.58 

-72. 5C 

192.84 






0 .42 196- C 1 

ie4.56 

-72.51 

194.30 






0.4265E-C1 

167.51 

-50. OC 

181.23 







the splitter blade. On the complete blade, there is a short distance after the end of the 
splitter blade before the complete blade assumes the blade loading without a splitter blade. 


Axial Stator 

This example is a stator nozzle mean blade section (fig. 7) for a turbine built at 
Lewis Research Center (ref. 9). This blade section has also been analyzed by using the 
incompressible flow program of reference 5. Downstream of the blade V/V cr is 0. 58. 
The input for this case is given in table II. The surface velocities obtained by the pro- 
gram are given in table III. The velocities are plotted against blade surface length in 
figure 8. Also shown in figure 8 are experimental data obtained from the investigation 
described in reference 9. There is fairly good agreement with the computed values for 
this example. 
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DESCRIPTION OF INPUT AND OUTPUT 


The computer program requires as input sufficient information to describe the blade 
shape accurately, the inlet and outlet angles, the extent of the region to be considered, 
the mesh size to be used, the appropriate gas constants, and operating conditions such as 
inlet temperature, density, weight flow, and rotational speed. Output obtained from the 
program includes velocity magnitude and direction at all interior points, blade surface 
velocities, stream function values, and streamline locations if there is no reverse flow. 

If there is reverse flow (as may occur with radial flow), streamline locations may be ob- 
tained by plotting contours of equal stream function values. 


Instructions for Preparing Input 

Figure 9 shows the input variables as they are to be punched on the data cards. 
There are two types of input variables, geometric and nongeometric. The geometric in- 
put variables are shown in figures 10 and 11. Units are always meters for length and 
radians for 9 coordinates. 



Figure 9. - Input form. (Card column numbers appear at top. ) 
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Direction of 
rotation 


o Spline point 



Figure 10. - Geometric variables required as input. Blade -to -blade coordinates on stream surface The 
variables BETAI, BETAO, ALUI, ALLI, ALUO, and ALLO are to be given as a true angle p, not the angle as 
measured in the m.Oplane. (Use tan p = rldB/dm) to obtain the value of p. ) 



Figure 11. - Geometric variables required as input. Meridional plane. 


The blade shape is defined by specifying the leading and trailing edge radii and a 
number of blade surface m and 9 coordinates. These coordinates are used to define a 
cubic spline curve (refs. 5 and 10). The coordinates are given with respect to the leading 
edge of the lower blade, as shown in figure 10. The standard sign convention is used for 
angles, as indicated in figure 10. The blade should be oriented with the concave side 
down. 

The mean stream surface of revolution and normal stream channel thickness are also 
defined by cubic spline curves as indicated in figure 11. The m coordinates for the 
mean stream surface are independent of the blade shape coordinates. 
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L 





HU 


Figured. - Mesh used for axial stator numerical example. Numbers are mesh point indexes (I in program). There are 840 unknown mesh points. 

A mesh must be used for the finite -difference solution of equation (1). A typical 
mesh pattern is shown in figure 12. The mesh spacing and the extent of the upstream and 
downstream regions are determined by the values of MXBI, MXBO, MX, and NBBI. The 
mesh spacing must be chosen so that there are not more than 2500 unknown mesh points. 

The values of /3. n and i3 Qut must be given on AH and DE, respectively. However, 
it may be that the average values along BG and CF are what is known. In this case the 
input values, /3. and P Qui , must be calculated by equation (B13) or (B15). 

While the program was written for compressible flow, it can be used equally well for 
incompressible flow. To use the program for incompressible flow specify GAM =1.5, 
AR = 1000, and TIP = 10 as input. Of course, RHOIP is simply the density in kilo- 
grams per cubic meter. This results in one iteration only. 

The International System of Units (ref. 11) is used throughout. However, the pro- 
gram does not use any constants which depend on the system of units being used (other 
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than degrees or radians where specified) so that any other consistent set of units may be 
used. For example, force, length, temperature, and time units may be chosen indepen- 
dently (mass units defined by f = ma); the gas constant R must then have the units of 
force times length divided by mass times temperature (energy per unit mass per degree 
temperature). Output then gives the velocity in the chosen units of length per unit time, 
and, of course, the output is not labeled with the correct units. 


Input 

All the numbers on the card beginning with MXBI and on the card beginning with 
BLDATA are integers (no decimal point) in a 5-column field (see fig. 9). Integers must 
be right adjusted. The remaining input variables are real numbers (punch decimal point) 
in a 10-column field. 

The input variables are as follows: 

GAM specific heat ratio, y 

AR gas constant, R, J/(kg)(°K) 

TIP inlet total temperature, T? n , °K 

RHOIP inlet total density, p| n , kg/m 3 

WTFL mass flow per blade for the stream channel, kg/sec 

OMEGA rotational speed, oo, rad/sec (Note that o> is negative if rotation is in the 

opposite direction of that shown in fig. 10.) 

W value for overrelaxation factor O to be used in eq. (A 8) (If W = 0, the pro- 

gram calculates an estimated value for the overrelaxation factor, see 
appendix A for discussion.) 

CHORD overall length of blade in m direction, meters, see fig. 10 

STGR angular coordinate 9 for center of trailing-edge radius, radians, see fig. 10 

BETAI flow angle j3- n along AH, deg, see fig. 10 

BETAO flow angle 0 out along DE, deg, see fig. 10 

RI leading-edge radius, m, see fig. 10 

ALUI angle at tangent point of leading-edge radius on upper surface, deg, see fig. 10 

ALLI angle at tangent point of leading-edge radius on lower surface, deg, see fig. 10 

RO trailing-edge radius, m, see fig. 10 

A LUO angle at tangent point of trailing-edge radius on upper surface, deg, see fig. 10 
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ALLO angle at tangent point of trailing-edge radius on lower surface, deg, see fig. 10 

MXBI number of mesh lines from AH to BG inclusive (fig. 12) 

MXBO number of mesh lines from AH to CF inclusive (fig. 12) 

MX total number of mesh lines in m direction from AH to DE, maximum of 100 

(fig. 12) 

NBBI number of mesh spaces in 0 direction between AB and HG, maximum of 50 

(fig. 12) 

NUSP number of blade spline points including end points that are tangent to leading- 

and trailing-edge radii for upper surface (BC) of blade (figs. 4 and 10); 
maximum of 50 

NLSP same as NUSP, but for the lower surface (GF) of the blade 

NRSP number of spline points for stream sheet radius (RMSP) and thickness (BESP) 

coordinates (see fig. 11), maximum of 50 

NBL number of blades 

NINT number of streamlines desired as output, maximum of 10 

MU array of m -coordinates of spline points for upper surface measured from 

leading edge, m, fig. 10 (The first and last points must be left blank, since 
these points are calculated by the program. If the last point is on a new 
card, a blank card must be used. The total number of points is NUSP.) 

XSPU array of 0-coordinates corresponding to the MU array, rad 

ML same as MU but for lower surface (The total number of points is NLSP.) 

XSPL array of 0-coordinates corresponding to the ML array, rad (Note that these 

coordinates are to the lower blade as shown in fig. 10.) 

MR array of m -coordinates of spline points for stream surface radii and stream 

channel thickness measured from leading edge, m, see fig. 11 (These 
coordinates should include the entire distance from AH to DE, and may ex- 
tend beyond these points if desired. The total number of points is NRSP.) 

RMSP array of stream surface radii corresponding to the MR array, m 

BESP array of stream channel thicknesses corresponding to the MR array, m 

The remaining variables, starting with BLDATA, are used to indicate what output is de- 
sired. A value of zero for any of these variables will cause the output associated with 
that variable to be omitted. A value of 1 will cause the corresponding output to be printed 
for the final iteration only; 2, for the first and final iteration; and 3, for all iterations. 
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Care should be used not to call for more output than is really useful. The following list 
gives the output associated with each of these variables. 


BLDATA 


NULAKI 


ERPRT 

STRFN 

SLCRD 

ARPRT 


INTVEL 

SURVEL 


radii and stream sheet thickness at each vertical mesh line; coordinates, 
first and second derivative of blade spline points; 9 -coordinate and slope at 
each vertical mesh line for each blade surface; coordinates of intersection 
of horizontal mesh lines with blade; NU and NL arrays (internal variables) 
(This will be printed for the first iteration only since these values do not 
change.) 

coefficient array A, the vector K, and the value of I for the adjacent points 
II, 12, 13, and 14 (This information is needed for debugging the program 
only.) 

the maximum change in the stream function for each iteration of the SOR 
equation, eq. (A 8) 

value of stream function at each unknown mesh point in the region 

streamline coordinates at each vertical mesh line and streamline plot 

values for (pW m ) and (pW 0 ) at all interior points and along blade surfaces, 
value of pW at all interior points (This information is needed for debugging 
the program only.) 

velocity and flow angle at all interior mesh points 

m -coordinate, surface velocity, flow angle, distance along surface, and pW 
based on meridional velocity components; m -coordinate, surface velocity, 
flow angle, and pW based on tangential components; plot of blade surface 
velocities against meridional streamline distance m (It is suggested that 
SURVEL = 3 be used. This will give surface velocities on every iteration, 
so that satisfactory velocities may be obtained even when final convergence 
is not reached in the alloted time.) 


Output 

Sample output is given for the radial inflow turbine numerical example, but with an 
outlet angle of -66. 5°. Since the complete output would be lengthy, only the first few 
lines of each type of output are reproduced here. Most of the output is optional and is 
controlled by the last input card as already described. The output labels are either in- 
ternal variable names or else are spelled out (e. g. , THETA for 9 ). 
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Each section of the sample output (table IV) has been numbered to correspond to the 
following descriptions: 

(1) The first output is a listing of the input data. All items are labeled as on the input 
form. 

(2) The calculated value of X is given followed by W (free -stream velocity) and the 
maximum value of the mass flow parameter pW (corresponding to W = W cr ) along AH 
(inlet) and DE (outlet). As a check, the free-stream values of 13 at leading edge (BG) and 
trailing edge (CF) corresponding to the input values of ( 3 ^ and /3 respectively, are 
calculated and printed out under the heading "BETA CORRECTED TO BLADE LE OR TE". 
The relative critical velocity W cr at BG (inlet) and CF (outlet) is printed out. Also 
some internal program constants are printed out at this point. 

(3) This is the output corresponding to BLDATA (see the list of input variables). 

(4) This is the number of mesh points at which the stream function is unknown. 

(5) This is the output corresponding to NULAKI. 

(6) If the program calculates an optimum overrelaxation factor £2 (i.e. , W = 0 for 
input) , then upper and lower bounds for £2 (WMAX and WMIN) and upper and lower bounds 
for p(Lj), (LMAX and LMIN) are printed out for each iteration (see appendix B of ref. 5 
for details). The last printed value of WMAX is the value of £2 (W) used by the program. 

(7) This is the output corresponding to ERPRT. 

(8) This is the output corresponding to STRFN. 

(9) This is the total execution time after obtaining the stream function solution for 
each outer iteration. 

(10) This is the output corresponding to SLCRD. 

(11) This is the output corresponding to ARPRT. 

(12) This is the output corresponding to INTVEL. 

(13) This gives the maximum relative change in the density p for each outer iteration. 

(14) This is the output corresponding to SURVEL. 

(15) This is the total execution time after all calculations are completed for an outer 
iteration. 
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Descrip- 

tion 

(a) 


TABLE IV. - SAMPLE OUTPUT 


GAM 

AR 

HP 

RH01P 

NT FL 

ONEGA 

U 


1.6667000 

208.20000 

1 C03.OOOO 

0.3956600 

0.62900 00 E-0 3 

*030.0000 

0 


CHORC 

STGR 

BETA1 

BET AC 





0 .684*000 E-C 1 - 

0.5350000 

-54.200000 

-66.500000 





RI 

ALU1 

ALLI 

RO 

ALUO 

ALLO 



0 .64 8000CE-03 

2. COCOOOO 

-2.0000000 

0 . 7 530 000 E-03 

-63.800000 

—63.800000 



MX 81 MX 80 MX NBB 1 

NUSP NLSP NR SP 

KBL NiNT 






5 39 *5 8 

11 11 16 

22 5 






MU ARRAY 








-0 

0. 86000 OOE- 02 

0.1600QGOE— 01 

0.2350000 E-01 

0.2903000 E-01 

0 .3*2 80 OOE- 01 

0. 395400QE— 01 

0. 4623000E-01 

0 .539*00 C E— 0 1 

C .6 1970C0E— 01 

-0 






XSPL ARRAY 








-0 

0.1220000E- 01 

0. 1 56U0C0E— 01 

0.1890000 £-01 

0.2090000 £-01 

0 .21 SOOOOE- 01 

0. 1600000E-01 

-0.11600 OOE-Ol 

-0 .96 20000E— 0 1 - 

C. 2791000 

-0. 542 7 COO 






ML ARRAY 








-0 

0. 06COOOOE— 02 

0. 1600000E-01 

0.2350000 £-01 

0.290 3000 E-01 

0 . 2A28CC0E- 01 

0. 3954000E— 01 

0.462 3000 E-01 

0.5394000^01 

C.6157000E-01 

-0 






XSPL ARRAY 







-0 

0.122000 OE— 01 

— 0. 1560000E-01 

-0.1900000 E-01 

-0.20 70000 E-0 1 

-0 . 222CCC0E- 01 

-0.2840000E-01 

-0.5 5 800 OOE-Ol 

-0 • 1* 5200 C 

0.33*7000 

-0. 5842 COO 






MR ARRAY 








—0 .762000 0 E-02 

0 

0. 8600000E-02 

0.1600000 E-01 

0.2350000 E-01 

0.25C3CCCE-01 

0.342 8000E-01 

0. 3 95 40 OOE -0 1 

0. *623000 E-01 

0.5394000E-01 

0. 61 97 0 OOE— 01 

0. 68440 OOE-Ol 

0.7354000E-01 

0 . 61 1 6CC0E- 01 

0.887 8000E- 01 

0. 9700000E-01 

RMSP ARRAY 







0.8A07000E— 01 

0 .76* 500 OE— 01 

0. 6800CC0E-O1 

0.6103000 E-01 

0.547 1000 E-C 1 

0 . 5G85CCCE- 01 

0.4808GOOE— 01 

0.46 02 000E-01 

0.4435000 E-01 

0.4295000E— 01 

0. 4131 G00E-01 

0 • hOO 5000 E-0 1 

0. 396 *000 E-0 1 

0 .3544GC0E-01 

0. 3 940000E — 01 

0.3 940000 E-01 

BESP ARRAY 







0. 9700000 E-03 

0 • 56 COO OOE— 03 

0. 1 03 0O00E-02 

0. 1090000E-02 

0.1140000 E-02 

0.116C0C0E-02 

0. 1 1600 OOE- 02 

0. 12400 OOE-02 

0.1330000 E-02 

0 • 1*3000 OE— 02 

0. 1530000E-02 

0.1620000 E-02 

0 .16 7 00 00 E-0 2 

0 • 16 5CCOOE- 02 

0.1 700000E-02 

0.1700000 E-02 

BLCATA NULAKI ERPR1 STRFN SLCRO ARPRT IMVEL 

SURVEL 





1 1 

3 2 

2 12 

3 






2 


L AHECA = 2C.136C0 

FREE STR EAM 
VELOCITY 

INLET 126.C0675 

OUTLET 250. £9178 


MAXIMUM VALUE 
FOR RHO * k 
125.45426 
105.286 52 


BETA CORRECTEC TO 
BLAOE LE OR TE 
-28. *27796 
-65.998662 


BLADE CRITICAL 
VELOCITY 
514.7898* 
497.7852* 


CALCULATED PROGRAM CONSTANTS 
HA HB HU 

0. 20 12 94 It— 02 C.356SS92E-01 C. 35012448-02^ 


HL 

0 • 3219 867 £-0 1 


PITCH 
0 .2855993 


N880 

-£ 


NBUO 


-16 


3 


STREAM SHEET COORDINATES AND THICKNESS TABLE 


M 


— 0 .805 1 8E— C2 
-0.60388E-C2 
-0.4025SE-C2 
-0 .20129E-C2 
0 


C . £* 5C4E- 01 
C.82482E-01 
C.8C466E-01 
C.7845AE-01 
C. 7645CE— 01 


SAL 
-1.00521 
-1. 0031 6 
-1.00060 
-0.99755 
-C. 95398 


0.97247E-03 
O.96235E-03 
0.9 56 3 IE— 0 3 
0.95523E— 03 
0.96000E-03 


CB/ DM 

—0 .58887 E— 02 
-0.40874E-02 
-0. 18435 E-02 
0.84291E-03 
0 .397 18E— 02 



Descrip- 

tion 

(a) 


TABLE IV. - Continued. SAMPLE OUTPUT 


3 


BLADE OATA AT SPLINE PCI MS 
UPPER SURFACE 

M THETA DERIVATIVE 

0.62754E-C3 C.E4719E-C2 C. -41324 

0 . 86000t- C2 C.122CUE-C1 C. 47782 

0 • 16000E- C 1 C • 1 5oCC£-0 1 C. 45426 


2ND DERIV. 

24.7234 

-8.52274 

2.15530 


3 


M 

0.62754E-C3 
0 .860GCE-C2 
C . 160o0t- C 1 


LUtaEK 
THETA 
C. 27713 
C. 27340 
C.2700C 


SURF ACE 

DERI VATI VE 
-o.41324 
-C. 47452 
-C.467C2 


2ND DERIV 
-25.5529 
10.1817 
-8.15426 


3 


BLAUE COORDINATE TABLE 


H 

0 

0.20129E-C2 

0.4G25SE-C2 

0.60386E-C2 

0.E0518E-C2 


XU 

0 

C. SC663E- 02 
C.95917E-C2 
C. 1C960E-01 
C.11937E-C1 


DAD2U 

50C0CCC0U000 
0.44349 
C. 47318 
0.48598 
0.48187 


XL 

0.2856U 
0.27653 
0.27560 
0.2 746 4 
0.27366 


3 


M 

0 

0.4230E-C 1 
0 .4424 E-C 1 
0.4814 E— C 1 
0 .4919E-C 1 


1NU 1NL 

1 


2 

3 


2 


1 TP 


1 

2 

2 

3 


4 { 


NUMBER CF INTERIOR MESH POINTS = 328 


3 


a See p. 20. 


CXUZL 

50UUO0000000 
-0 .44434 
-0 .47420 
-0.48589 
-0.47943 
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TABLE IV. - Continued. 


SAMPLE OUTPUT 


Descrip- 

tion 

(a) 


M H AX 

* 

2.O00CCO 

MM IN 

« 

1 . OCC 653 

LMA X 

a 

1.000000 

LM1N 

a 

0 .002768 

MMAX 

* 

2 . COCOCO 

KM IN 

* 

1.125C51 

LMAX 

a 

1.000000 

LM IN 

a 

0 .395284 

MM AX 

* 

2.UCCCCC 

taMIN 

= 

1.2CC5E6 

LMAX 

= 

1.000000 

LMIN 

a 

0 .557378 

WM AX 

M 

1.99C551 

MM IN 

a 

1.234761 

LMA X 

a 

0.999980 

LM1N 

= 

0.615915 

UMAX 

* 

1.56C779 

MM IN 

= 

1.246855 

LMAX 

a 

0.999600 

LMIN 

a 

0 .635141 

MM AX 

= 

1.950 149 

MH IN 

a 

1.250514 

LMA X 

= 

0.999347 

LMIN 

= 

0.640789 

MM AX 

= 

1.535527 

MM I N 

= 

1.261752 

LMAX 

= 

0.999041 

LMIN 

= 

0.657660 

.. UMAX 

* 

1.930263 

MM IN 

» 

1.274517 

LMAX 

3 

0.998695 

LMIN 

- 

0.675987 


7 


ERROR * 1.03 <56 773 
ERROR = 0.66 680200 
ERROR » 0.8EE3522S 
ERROR - 0. 41*45254 
ERROR = 0.56563662 
ERROR =* 0.25254514 
ERROR = C .38 5245C9 


- 

STREAM 

FLNCTIUN VALUES 






IA 

= 1 

-C. 35951178 

-C. 23676675 

-0.11125914 

0.01617927 

0.14364282 

0.2697986C 

C.3543C712 

0.51764261 

IA 

3 2 

-C. 243866*6 

-0.12112157 

C.C0436607 

0.13182450 

0.25928804 

0 .38544382 

C.5C555234 

0.633287 82 

IA 

3 3 

-C. 14629163 

-0.C2469553 

0.1C10C400 

G. 22934018 

0.35743163 

0.48374735 

C.6C603247 

C. 73104610 

IA 

= 4 

— C .065656 1 7 

G .C 5 125 642 

C. 17862408 

0.30920815 

0.43870088 

0.5653C646 

C. 66856664 

0.81121481 

I A 

C. 104443*4 

C. 2366*545 

C. 3 72C2 834 

0.50390199 

0 • O3082559 

0.75311223 

0.6 734C546 
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{ TIME = 0.2e26 MIN. 


STREAMLINE COORDINATES 


10 


M COURC. 

STREAM FN. 

THE TA 

STREAM FN. 

TFETA 

STREAM FN. 

THETA 

-C.8C51765E-C2 

— C ■ 20GC0CC 
C.4000COO 

C. 4624002 E—0 1 
0.2158414 

0 

G.6G0U000 

0 .1C25766 
0.2738242 

0.2C00000 

0.1586906 

— C.6C38E23E— C2 

— C . 20CCOCG 
C . 4 GCCC uO 

C. 12 83 459 E— 0 1 
0.1826500 

0 

0.6000000 

C .?C 16416 E— Cl 
G. 2402432 

C . 2CC0000 

0.1261657 

-C.4025862E-C2 

-c 

C • 6000CCO 

C. 42 858 18 E -0 1 
0.2118780 

0 .2000000 
0.8000000 

G.S6955C7E— Cl 
0.2655632 

0.40C0000 

0. 1547645 

-0.2U12S*1E-C2 

-c 

C. 60COCUO 

G.2040422E-0I 

0.1884405 

0.2000000 

0.8000000 

0.7 725C71E-C 1 
0.2466263 

U.40C0000 

0.1320617 

0 

0 

0 

U.200U000 

C. 6155541 E— Cl 

C.4C00000 

0.1146033 


C . 600GG0O 

C. 1696966 

O.BOOOUOO 

0 .2281752 

1.0000000 

0.2 855993 


a See p. 20. 
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TABli: IV. - Continued. SAMPLE OUTPUT 


STREAMLINE FLCTS 


0.042 

1 * — 


0.125 

* — 1 * — 

1 


0.291 

*-l — 

1 


STREAMLINES ARE PLOTTED WITH THETA ACROSS TEE PAGE AND M DO MN THE PAGE 


a See p. 20. 
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Descrip- 

tion 

(a) 


TABLE IV. - Continued. SAMPLE OUTPUT 


VELOCITIES AT INTERIOR MESH POINTS 




I A* 1 

VELOCITY 

ANGLE (DEG) 

VELOCI TY 

ANGLE (CEO) 

VELGC1TY 




131. 75 

-57.02 

133.12 

-56.51 

133.96 




132.66 

— 56.CS 

132.41 

-56.43 

132.25 

12 

4 

I A* 2 

VEL JCIT Y 

ANGLE (DEG) 

VELOCI TY 

ANGLE (DEG) 

VELOCITY 




121.41 

-52.11 

122.18 

-51.33 

123.25 




1 23.25 

-51. 3C 

122.57 

-51.61 

122.16 


ANGLE1 DEG) 

VELOCITY 

ANGLE 10EG) 

VELOCITY 

ANGLE (DEG) 

-55.96 

-56.63 

134.01 

-55.66 

133.49 

-55.76 

ANGLE! OEGJ 

VELOCITY 

ANGLE <OEG) 

VELOCI TY 

ANGLE (DEG) 

-50.78 

-51.75 

123.56 

-50.65 

123. 85 

-50.91 


I A ^ 


V. 


VELOCITY 

ANGLE (OEG ) 

VELOCI TY 

ANGLE IDEG) 

VELOCITY 

ANGLE! OEG) 

VELOCITY 

ANGLE (DEG) VELOCITY 

ANGLE (OEG) 

11U.54 

-46.77 

111.07 

-44 . 7b 

113.49 

-44.26 

114.55 

-44.52 115.02 

-45.11 

1 14.2C 

-45.66 

113.13 

-45.95 

112.42 

-45.94 





13 { ITERATION NO . 


MAXIMUM RtLATlVE CHANGE IN DENSITY = 0.8677 


14 





* 


SURFACE VELOCITIES 

BAS EC ON 

AX (AL 

COMPONENTS 






* 


UPPER SURFACE 


+ 


LOWER SURFACE 



M 

* 

VELOCITY 

ANGLE (DEG) 

SLRF. LENGTH 

KHC*W 

• 

VELOCITY 

ANGLE (OEG) 

SURF. LENGTH 

RHG*W 


0 

♦ 

0 

SC. oc 

0 

20.954 

• 

0 

-50. CO 

0 

31.393 


0 .2G13E-C* 

* 

12C.56 

1.85 

0.2126E-02 

42.299 

+ 

75.725 

-1.56 

C. 21266- 02 

28.303 


0 .4026E-C2 

* 

138.58 

1.96 

C.414GE— 02 

4B.095 

» 

57.442 

-2.C4 

0.4140E-G2 

20.343 


0 .6C3SE-U2 

* 

isi.es 

1.96 

0.6154E— 02 

51.912 

+ 

46.585 

-2.C3 

0.6154E-02 

16.418 


0 .EU5*E-C2 

* 

161. C4 

1. 85 

0. 81680-02 

54.465 

+ 

41.455 

- 1.55 

0.8168E-02 

14.539 


0 . 1006 £— C 1 

* 

16 7.66 

1.78 

C.1U18E-01 

56.167 

+ 

35.853 

-1.E2 

0. 101 BE — 01 

13.894 



SURFACfc V EEOC 1 I IE S BASED ON TANGENTIAL 
UPPER SURFACE 

M VELOCITY ANGLE IUEGI 


COMPONENT S 
RHC*h 


14 


+ 02 * 


VALUE OF RHU*M IS TOO LARGES 
CICK NAHe IFN OF CALL 
.ERROR CC025 

CENSIY QUO 1 6 

TASVEL OG 524 

2 COP GOO 1 8 


ABS. LOC. 
16615 
1C132 
2ol 1 3 
03056 


-0 53.341 

0.4424E-C1 184 .SC 

0.4919E— C 1 19 5.52 

0.5227E-C 1 206. 3 fc 

0 • 5460 fc— C 1 2 1 £ . 32 

0 .56 5 2E- C i 231. 5C 


90. GC 19.164 

-12.61 57.795 

-23.06 6C.602 

-30.50 63.453 

-36.16 66.525 

-40.57 69.822 


14 


M 

0.423UE-C 1 
0.4814E-C 1 
C.5133E-C 1 
0 .5379E-C 1 
0.5584 E-C 1 
„ 0.5759E-C 1 


VELOCI 
129.95 
11C. 54 
1 19 • £5 
12 £ . 55 
135.65 
143.75 


LONER SURFACE 
TY ANGLE (DEG 1 
-8. 87 
-22. 23 
-25. 54 
-34. 52 
-38. 86 
-43.13 


RHC *h 
41.775 
35.681 
38. 524 
41.130 
43.255 
45.644 


a See p. 20. 


25 



i 


Descrip- 

tion 

(a) 


i \ 


TABLE IV. - Continued. SAMPLE OUTPUT 


BLAOE SURFACE VELOCITIES 


0. 

so. 

LOO. 

150. 


200. 


25C. 

3CC. 

350. 

40C. 

450. 

500. 

1 

X 

♦ X 

* 1 


1 



1 

1 

1 

1 

1 

l 

1* 

X 

* 1 


1 


1 

1 

1 

1 

1 

1 

1 


X 

* 


1 


1 

1 

1 

1 

1 

1 

I 

♦ 1 

1 

1 * 


1 


1 

1 

1 

1 

1 

1 

1 

♦ X 

1 

1 

• 

1 


I 

1 

1 

1 

1 

1 

X 

♦ X 

1 

1 

* 

1 


1 

1 

1 

1 

1 

1 

1 

+x 

1 

1 

* 

1 


1 

1 

1 

1 

1 

1 

1 


L 

1 

* 

1 


1 

1 

1 

1 

l 

1 

X 

X* 

1 

1 

* 

i 


1 

1 

1 

1 

1 

1 

C.C2CX 

— . — !-♦ — 

X 

1 — 

— 

— 1 

— 

— 1 

— 1 

1 

1 

1 — 

1 — ■ 

X 

X ♦ 

1 

1 

♦ 

1 


J 

1 

1 

1 

1 

1 

1 

X 

«• 1 

1 

♦ 

1 


1 

1 

1 

1 

1 

1 

X 

X 

+ 1 

1 

* 

1 


1 

1 

1 

l 

1 

1 

X 

X 

+ 1 

1 

* 

1 


1 

1 

1 

1 

1 

1 

X 

X 

♦ 1 

1 

* 

1 


1 

1 

1 

1 

1 

1 

X 

1 

1 + 

1 

• 

1 


1 

1 

1 

1 

1 

1 

X 

X 

1 ♦ 

1 

♦ 

1 


1 

1 

1 

1 

1 

1 

X 

1 

1 + 

1 

* 

1 


1 

1 

1 

1 

1 

1 

X 

X 

X * 

1 

* 

1 


1 

1 

1 

1 

1 

1 














1 

1 

I 

1 

• 

1 


1 

1 

1 

1 

1 

1 

1 

X 

X* 

1 

* 

1 


1 

1 

1 

1 

1 

1 

1 

X 

1 * 

1 


* 1 


1 

1 

1 

1 

1 

1 

X 

X 

X «■ 

1 


♦1 


1 

1 

1 

1 

1 

1 

X 

X 

1 ♦ 

1 


* 


1 

1 

1 

1 

1 

1 

X 

X 

1 

♦ 1 


1 * 


1 

1 

1 

1 

1 

1 

1 

X 

1 

♦ 1 


1 *0 


1 

1 

1 

1 

1 

1 

1 

X 

1 

♦ 1 


1 

0 

1 

1 

1 

1 

1 

1 

1 

X 

1 

X* 


1 


C* 

1 

1 

1 

1 

1 

X 

X 

I 

1 

X 

♦XI 


1 

* 

1 

1 

1 

1 

I 

X 

X 

1 


1 X 

X 

1 

1* 

1 

1 

1 

1 

l 

1 

1 

1 


1 


X X 

1 00 c 

) 1 

1 X 

1 

1 

X 

1 

1 

1 


1 


1 

1 C 

01 

1 

1 

0 

1 

l 

X 

1 


1 


1 

1 

1 

1 

1 

1 

X 

1 

1 

1 


1 


1 

1 

1 

1 

1 

1 

X 

X 

1 

1 


1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 


i 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 


1 

1 

1 

1 

1 

1 














c. 

50. 

100. 

150. 


200. 


25 C. 

3CC. 

350. 

400. 

450. 

500. 


VELCCI TY IMETERS/SECLNG) VS . HER 10 IONAL STREAMLINE C I S TANCE { ME Tfc R S ) DOWN THE PAGE 


0 

X 


UPPER SURFACE, 
LOWER SURFACE, 
UPPER SURFACE, 
LGWER SURFACE, 


BAS EC CN 
BASED ON 
BASED ON 
BASED CN 


AXIAL COMPONENT 
AXIAL COMPONENT 
TANGENTIAL COMPONENT 
TANGENTIAL COMPONENT 


^ee p. 20. 


26 



Descrip- 

tion 

(a) 

15 { T IME - G .35C6 MIN. 


TABLE IV. - Continued. SAMPLE OUTPUT 


7 


ERROR ■ 
ERROR » 
ERROR - 
, ERROR - 
ERROR - 
ERROR - 
ERROR * 
„ ERROR « 


0 . 0 1299 141 
C.CC540415 
0 .GC 155128 
0 .OC 6C3 5C 1 
0.0C576552 
0.00 565254 
0.0C £C53C3 
C.OC ‘38156 


9 ^ T IME - 0 .56 CC MIN . 


{ ITERATION NO. 2 MAXIMUM RELATIVE CHANGE IN DENSITY * 0.3426E-01 


SURFACE VELLCITltS 8ASEU CN AXIAL LUMRONENIS 
SURFACE 


M 

* 

VELOCITY 

ANGLE (OEG) 

SLRF. LENGTH 

RNG*W 

• 

VELOCITY 

ANGL E (DEG ) 

SURF. LENGTH 

RHG* W 

0 

* 

0 

50. OC 

0 

16.887 

* 

0 

-50.00 

0 

37.242 

0 . 20 13E-0 2 

* 

1C6.83 

1.85 

0. 212 6E-02 

37.651 

* 

94.635 

-1.56 

0.2126E-02 

33.541 

0 .4U26E-0 2 

* 

125.61 

1.96 

0 .41406— 02 

43.692 

• 

69.861 

-2.C4 

0.4140E-02 

24.665 

0 .6039E-G2 

* 

126.13 

1.96 

0 . 6154E-02 

47.498 

• 

58.054 

-2.C3 

C.6154E-02 

20.423 

0.8052E-C2 

• 

146. 8C 

1.85 

0. 61o8E-02 

49.972 

* 

52.6<6 

-1.55 

0 . 81 686 — 02 

16.412 


SURFACE VELOCITIES BASEO ON TANGENTIAL 
UPPER SURFACE 

M VELOCITY ANGLE (DEG) 


COMPONENTS 

RHC*W 


14 


* 04 * 




ABS. LCC. 
16615 


0ENS1Y GOO 1 fc 

TASVEL 00524 

2LCP 0CC18 


1C132 

26113 

C3C56 


•0 12.175 

0.4424E-C) 176.45 

0.4919E-C1 165.28 

0.5227E-C1 20C.12 


50. CC 25. 845 

-12.61 55.993 

-23. 06 58.913 

-30. 5C 61.772 


14 


M 

0 .42306— C 1 
0 .48146— C 1 
0 . 5 1 3 3 E— 0 1 
0.5375E-C 1 


LUMER SURFACE 
VELOCITY ANGLE (DEG) 
135. 7E -8.87 

114. *5 -22.23 


123.55 -25. 54 

132.23 -34. 52 


PHC*t* 

43.533 

36.842 

35.647 

42.246 


‘See p. 20. 



TABLE IV. - Continued. SAMPLE OUTPUT 


Descrip- 

tion 

(a) 


14< 


1* 

1 + 
1 

0 * 0201 — 
1 
1 
X 
1 
1 
1 
1 
1 
1 

O.C401 

1 

1 

1 

1 

1 

1 

1 

1 

1 

C.C6CL 

1 

1 

1 

1 

1 


BLADE SURFACE VELOCITIES 
. 200 . 250 . 

i 1 

1 1 

1 1 


1 

1 

1 

- 1 — 

1 

1 

1 

♦ 

1 + 

1 ♦ 
1 
1 
1 

X ♦ 
1 * 
1 + 
1 
1 
1 
1 
1 
1 

1 

1 

1 


i* 

1 * 
1 
1 

-1 — 
1 
1 
1 


1 

1 

1 

• 1 

♦ 

1 i 

I 

1 

— 1 - 
♦ X 
1 
I 


I 

1 

1 

1 

1 

l 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


-l-o- 

1 

1 

IX 

1 

1 

1 

1 


I 

1 

1 

— 1- 
1 
1 
1 
1 
1 
1 

1 

I 

<1 
♦ 1 

cc 

c 

1 

1 

1 

1 

1 


450. 
— i — 


45C. 


500. 

1 


--I — 
500. 


VELCCI TV CPETERS/SECCNGI VS. HER ID 1 JNAL STREAMLINE C I S TANCE I ME IE R S > OChN THE PAGE 


550. 

1 


- UPPER SURFACE, 

- LGHER SURFACE, 

- UPPER SURFACE, 

- LCwER SURFACE, 


eASEC CN AXIAL COMPONENT 
BASED GN AXIAL CUMPUNENT 
BASEC CN TANGENTIAL CUMPUNENT 
BAS EC GN TANGENTIAL CUMPUNENT 


l See p. 20. 



TABLE IV. - Continued. SAMPLE OUTPUT 



Descrip- 

tion 

(a) 

15 { TIME - 0.64SE MIN. 




• 


SURFACE VELCCITIES 

BASED ON 

AX 1AL 

COMPONENTS 






* 


UPPER SURFACE 


♦ 


LOWER SURFACE 



M 

* 

VELOCITY 

ANGLE I0EGI SURF. LENGTH 

RHG*W 

* 

VELOCITY 

ANGLE IOEG I 

SURF. LENGTH 

RHO+W 


0 

• 

0 

50.00 0 

16.975 

* 

0 

-90.00 

0 

37.167 


0 .2013E— C2 

* 

107.31 

1.85 0. 2126E-02 

37.814 

* 

94.57E 

-1.56 

0.2126E-02 

33.451 


0.4026E-02 

* 

126.22 

1.96 0.4140E-0 2 

43.893 

* 

69.651 

-2.C4 

0.4140E-02 

24.612 


0.6039E-02 

* 

138.86 

1.96 0. 6154E-02 

47.734 

* 

57.862 

-2.C3 

0.6154E-C2 

20.356 


L 0.e052E-C2 

* 

147.64 

1.85 0.8168 E— 02 

50.240 

* 

52.435 

-1.55 

0.8168E-02 

18.347 




14 


LOWER SURFACE 

H VELOCITY ANGLE (OEGI PHC*W 

0 . 4230E- C 1 135. At -0.87 43.436 

0.4814E-C1 114. C7 -22.23 36.775 

0.5133E-C 1 123.32 -25.54 39.573 

0.S379E-C1 131.56 -34.52 42.166 

0 .S5B4E-C 1 135. 2C -38. 86 44.305 


15 { TINE - C.85C6 MIN. 

^ee p. 20. 


29 


TABLE IV. - Continued. SAMPLE OUTPUT 


Descrip- 

tion 

(a) 

f BLADE SURFACE VELOCITIES 


SC. 

100. 

ISO. 


200 . 

250. 

3C0 . 

3 EC . 

4CG. 

45 C. 

500. 

550. 

1 

♦ 1* 

i 


1 

1 

1 

1 

1 

1 

1 

1 

l ♦ 

1 

♦ 1 


1 

1 

1 

1 

1 

1 

1 

1 

1 ♦ 

1 

♦ 1 


1 

1 

1 

1 

1 

1 

1 

1 

* 

1 

* 


1 

1 

1 

1 

1 

l 

1 

1 


1 

1* 


1 

I 

1 

1 

1 

1 

1 

1 

+ 

I 

1 • 


1 

1 

1 

1 

I 

1 

1 

1 

1+ 

1 

i * 


1 

1 

1 

1 

1 

1 

1 

1 

l + 

1 

1 

* 

1 

1 

1 

1 

1 

1 

1 

1 

1 ♦ 

1 

1 

+ 

1 

1 

1 

1 

1 

1 

1 

1 

1 * 

i 

1 

* 

1 

1 

1 

1 

I 

1 

I 

1 

1 

* I 

1 

• 

I 

1 

1 

1 

1 

1 

1 

1 

l 

♦ 1 

1 

* 

1 

1 

1 

1 

1 

1 

1 

1 

i 


1 

* 

1 

1 

1 

1 

1 

I 

1 

1 

1 

1* 

1 

* 

I 

1 

1 

1 

1 

1 

1 

1 

1 

1 ♦ 

1 

* 

1 

1 

1 

1 

1 

1 

1 

1 

I 

1 ♦ 

1 

* 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 ♦ 

I 

♦ 

1 

I 

1 

1 

1 

1 

1 

1 

1 

1 ♦ 

1 

♦ 

1 

1 

1 

1 

l 

1 

I 

1 

1 

1 ♦ 

1 

• 

1 

1 

1 

1 

1 

1 

1 

1 

I 

1 ♦ 

1 

• 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 + 

1 

* 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 ♦ 

1 

* 

1 

1 

1 

1 

l 

1 

1 

1 

1 

1 ♦ 

1 


♦1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

♦ I 


♦ 

1 

1 

1 

1 

1 

1 

1 

1 

1 

X+ I 


1 *0 

1 

1 

1 

1 

1 

1 

1 

1 

1 

+ 1 


1 

0 1 

I 

1 

1 

1 

l 

1 

1 

1 

X* 


1 

0*1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

X 

XI 

1 

* 

1 

1 

1 

1 

I 

1 

1 

1 


1 X 

X 1 

* 

1 

1 

1 

1 

I 

1 

I 

1 


I 

IX 

X OOC 

1 

1 X 

1 

1 

1 

1 

1 

1 


1 

1 

0 1 C 

1 

1 

1 

0 

1 

1 

1 

I 


1 

1 

1 

1 

1 

1 

1 

1 

I 

1 

1 


i 

1 

1 

1 

1 

1 

1 

1 

1 

I 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

I 

1 

i 


1 

1 

1 

1 

1 

1 

1 

1 

5C. 

100. 

150. 


200. 

250. 

3C0. 

3 EC • 

4CC. 

450. 

500. 

550. 


VELOCI T V ( METERS/S ECCNO ) VS. MERIDIONAL STREAMLINE D ISTANCE ( ME TER S ) DORN THE PAGE 


0 

X 


UFPER SURFACE, 
LOWER SURFACE, 
UPPER SURFACE, 
LOWER SURFACE, 


BASED CN AXIAL COMPONENT 
BASED ON AXIAL COMPONENT 
BASED CN TANGENTIAL COMPONENT 
BASED CN TANGENTIAL COMPONENT 


l See p. 20. 


TABLE IV. - Continued. SAMPLE OUTPUT 


Descrip- 

tion 

(a) 


5 


IA 

I 

A( 1. 1) 

All. 2) 

AU .3) 

All ,4! 

11 

12 

13 

14 

K ( 1) 

1 

1 

C. 

0. 

0. 

l .oooou 

8 

2 

0 

5 

-C. 11565 

1 

2 

c. 

C. 

c. 

1.00000 

1 

3 

0 

1C 

-C. 11565 

1 

3 

c. 

C. 

0. 

1.00000 

2 

4 

0 

11 

-C .1 1565 

1 

4 

c. 

C. 

c. 

l.OOOUO 

3 

5 

0 

12 

-C.l 1565 

1 

5 

c. 

c. 

c. 

1.00000 

4 

6 

0 

13 

-0.11565 

1 

6 

c. 

G. 

c. 

1.00000 

5 

7 

0 

14 

-C. 11565 

1 

7 

c. 

0. 

c. 

1.00000 

6 

8 

c 

15 

-C. 11565 

1 

E 

c. 

0. 

c. 

1.00000 

7 

1 

0 

16 

-C.l 1565 

2 

9 

C. 15527 

C. 15521 

0.34271 

0.33881 

16 

10 

1 

17 

-C. 15310 

2 

10 

C. 15524 

C. 15524 

C. 342 78 

C. 33874 

9 

11 

2 

18 

C.CC617 

2 

11 

C. 15522 

C. 15526 

0.34278 

C. :»38 74 

10 

12 

3 

15 

C.CC617 

2 

12 

C. 15522 

C. 15526 

0.34275 

C. 33877 

11 

13 

4 

2C 

C.GC616 

2 

13 

C. 15923 

C. 15525 

0.342 73 

0.33879 

12 

14 

5 

21 

C.CC616 

2 

14 

C. 15924 

C . 15524 

C. 342 72 

0.33880 

13 

15 

6 

22 

O.OC616 

2 

15 

C. 15524 

C. 15924 

C. 342 72 

G. 33880 

14 

16 

7 

23 

C .CC 6 16 

2 

16 

0.15527 

C. 15522 

C. 34269 

0. 33bs3 

15 

9 

8 

24 

C. 16 538 

3 

17 

C. 16472 

C.1646C 

C. 33815 

C. 33252 

24 

18 

9 

25 

-U.15fc73 

3 

18 

C. 16466 

G. 16466 

0.33820 

0.33248 

17 

19 

10 

26 

C .CC5S5 

3 

15 

C. 16462 

C.1647G 

C. 33815 

0.33253 

18 

20 

11 

27 

C.CC555 

3 

*0 

C. 16463 

C. 16469 

C. 33811 

C. 33257 

19 

21 

12 

28 

G.CC555 


” ERROR * 
ERROR « 
ERROR = 
ERROR = 
ERROR * 
ERROR = 
„ ERROR = 


O.uCCC.6517 
O.OCCC3825 
C.OCCC3228 
0.CCCC3EE5 
G.UCCC2SC3 
C .OC CC2 3 64 
C.uC CC 12 64 


STREAM FUNCTION VALUF5 


IA 

* 1 

-C .Jd5uuE45 

-G .263 7G5fc2 

-0.14113821 

-0.01622819 

0.11014967 

0.23700804 

C.3637S162 

0.49012671 

IA 

* 2 

-C. 26936 323 

-0 .14E0606C 

-0. C2545295 

C. 09941702 

0.22579490 

0.35265326 

0.47542664 

0.6C577193 

IA 

= 3 

— C • 16 869 7 2 C 

-0.C4983827 

C.C 7221C28 

C. 19751996 

0.32436176 

0.45148222 

C.57fc4<l35 

C.7C525534 

IA 

= 4 

-C. 08133*51 

C.C3C5756C 

0. 15i5CE31 

C. 27847333 

0.40642122 

0.53403744 

C.66 1 1 6 * 4 £ 

G.7888E537 

I A 

= 5 

C.C9U1U2 31 

C .2 1272253 

C. 3426CCU 

0.47260920 

0 .buO 83592 

0.72757553 

C .fc 5 54 7 2 1 £ 



9 { TIME = 2. 7744 M IN . 


STREAMLI NE 


C00KU1 NATES 


M LOORC. 


STREAM FN. 


IHETA 


ST REAM 


FN. 


TFETA 


STREAM FN. 


THETA 


-C.dOSi '/<££- C 2 
-C . 60 3d £ 23E- C2 
-C.'tC^bet 2E-C2 
-C.2Q12541E-C2 


-C. 2CC0M.C 
G. 4CC000G 
-C.2CLGGCC 
C. 4GC00C0 

-c 

C. 6COUC00 

-c 

O.bOCLOOC 

c 


C. 543269UE -0 1 

0.224<*047 

G.2044267E-C1 

G. 1*918262 

C.5CJ9991E-G1 

C.22U2t>U8 

C.2620723E-01 

C.197C188 

C 


0 

0 .6000000 
0 

0.6U00O00 
0 . 2U0GG00 
J . 6 JOUUOO 
U . 20C00CU 
0 . 8G0000G 
0.2000000 


0 . 1 ) 16555 
0.2612534 
G.7£?2S75E-C1 
0 .2482602 
C . 1C ?eCC2 
G. 2767515 
C.E5C6C24fc— C 1 
0 • 252S6C7 
C.tise 16 fc E- C 1 


G • 2CCOOOC 
0.2CC0000 
C.4CCC000 
0 . 4CCGC0C 
0.40C0Q0Q 


0. 1680872 
0.1355313 
0. 1640383 
0.1410080 
0.1228669 


a See p. 20. 
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TABLE IV. - Continued. SAMPLE OUTPUT 

Descrip- 

tion 

(a) 


r STREAMLINE PICTS 


-0.539 


-C.456 

-0.373 

-0.290 

-0.207 

-0.124 

-C.C41 

0.042 

0.125 

0.208 

0.291 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

* 1 * 

1* 

* 1 

* l 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

1* 

♦ 1 

* 1* 

* 1 

i 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

* 1 

1 * 

*1 

* 1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

* 1 * 

* 

* 1 * 

♦ 1 

1 


1 

l 

1 

1 

1 

1 

1 

1 

1 

1 

C. 00201 — 

— 

1 

1- — 

-1—— 

1— 

— 1- 

— 1 — -- 

— * 1-* 

— *-l 

* 1-*— 

*-l — 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


l 

1 

1 

1 

1 

1 

* l* 

* 1 * 

* 

* 1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

* 1* 

* 1 * 

*1 

* 1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

♦ 1* 

* 1 * 

*1 

* 1 

1 


1 

1 

1 

1 

1 

1 

1 

l 

1 

1 

1 


1 

1 

1 

1 

1 

1 

* 1* 

♦ 1 * 

*1 

* 1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

0.01201 

— 

1 

1 

1 

— - — 1 

1 

1 

* — 1* 

* 1-*— 

*-l 

*-l — 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

♦ 1* 

* 1 * 

* 1 

* 1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

* 1* 

* 1 * 

* 1 

♦ 1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

* 1* 

* 1 * 

*1 

* 1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

l 

1 


1 

1 

1 

1 

1 

1 

♦ 1 * 

♦ 1 * 

♦1 

* 1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

1 

l 

1 

1 

1 


1 

1 

1 

l 

1 

1 

* 1 ♦ 

* 1 ♦ 

*1 

* 1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

l 

1 

1 

1 

1 

* 1 * 

* l * 

*1 

* 1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

i 

1 

* 1 * 

* 1 * 

* 

* 1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

* 1 * 

* 1 * 

* 

* 1 

1 


1 

1 

l 

1 

1 

1 

1 

1 

1 

1 













1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

* l » 

* 1 * 

• 

* 1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

i 


1 

1 

l 

1 

1 

1 

* 1 * 

♦ 1 * 

* 

* 1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

* 1 * 

* 1 » 

*1 

* 1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

* 1 * 

* 1 * 

*1 

* 1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


I 

1 

1 

i 

1 

l 

1 

1 

i 

1 

1 


1 

1 

1 

l 

1 

1 

* * * 

1* 

* i * 

1 

1 


1 

1 

1 

l 

1 

1 

1 

1 

l 

1 

1 


1 

1 

1 

l 

1 

1 * 

* 1 * 

*1 

* i * 

1 

1 


1 

1 

1 

i 

1 

1 

1 

1 

l 

1 

1 


1 

1 

1 

l 

1 

1 * 

* 1 * 

* 1 * 

i* 

1 

1 


1 

1 

1 

i 

1 

i 

1 

1 

i 

1 

1 


1 

1 

1 

l 

1 

* l 

* *i * 

1* 

* i 

1 

i 


1 

1 

1 

l 

1 

l 

1 

1 

i 

1 

0.05201 

— 

— 1 

1 

1 

1 — 

1 

-* — l* — 

— * 1-* 

— *— 1 

* i 

— — — 1 — 

1 


1 

1 

1 

l 

1 

l 

1 

1 

i 

1 

1 


1 

1 

1 

l 

1 * 

* i * 

* 1 * 

1* 

i 

1 

1 


1 

1 

1 

l 

1 

l 

1 

1 

i 

1 

1 


1 

1 

1 

l 

* 1 * 

* l 

* * 

* 1 

i 

1 

1 


1 

1 

1 

l 

1 

i 

1 

1 

i 

1 

1 


1 

1 

1 

i * 

* 1 * 

* l 

* i* 

1 

i 

1 

1 


1 

1 

1 

l 

1 

l 

1 

1 

i 

1 

1 


1 

1 

1 

♦ l * 

* 1 * 

* l 

* 1 

1 

l 

1 

1 


1 

1 

1 

l 

1 

l 

1 

1 

i 

1 

1 


1 

1 

1 

l 

1 

} 

1 

1 

i 

1 

1 


1 

1 * 

* 1* 

* l • 

* 1 

l 

1 

1 

i 

1 

1 


1 

1 

1 

l 

1 

l 

1 

1 

i 

1 

1 


i * 

* 1 ♦ 

* 1 * 

* l 

1 

i 

1 

1 

i 

1 

1 


1 

1 

1 

i 

1 

l 

1 

1 

i 

1 

* 

* 

1 * 

• 1 

* 1 * 

l 

1 

l 

1 

1 

i 

1 

1 


1 

l 

1 

l 

1 

l 

1 

1 

l 

1 

1 

* 

1 * 

*1 

* 1 * 

i 

1 

l 

1 

1 

i 

1 

1 


1 

1 

1 

l 

1 

l 

1 

1 

l 

1 

0.07201 

— *- 

— 1-* 

— *1 

— * 1--*- 

1 — 

1 

1 — 

1 

1 

i 

1 — 

1 


1 

1 

1 

l 

1 

i 

1 

1 

i 

1 

1 

* 

l • 

*1 

* 1 * 

l 

1 

i 

1 

1 

i 

1 

1 


1 

1 

1 

l 

1 

i 

1 

1 

i 

1 

1 

* 

1 * 

*1 

* 1 * 

i 

1 

1 

1 

1 

i 

1 

1 


1 

1 

1 

l 

1 

1 

1 

1 

i 

1 

1 


1 

1 

1 

i 

1 

1 

1 

1 

i 

1 

1 

* 

1* 

* 1 

* 1 * 

l 

1 

1 

1 

1 

i 

1 

1 


1 

1 

1 

l 

1 

1 

1 

1 

i 

1 

1 

• 

1* 

* 1 

* 1 * 

l 

1 

1 

1 

1 

l 

1 

0.08201 

— 

— 1 

1 

1 

1 — 

1 

1- 

1 

1 — 

1 _ — 

1 — 

-C.535 


C.456 

-0.373 

-0.290 

-0.207 

-C.12< 

-C.C41 

C.042 

0.125 

0.208 

0.291 


STREAMLINES ARE FLCTTEC WIFh THETA ACROSS TEE PAGE AND N DOWN THE PAGE 


l See p. 20, 


Descrip- 

tion 



TABLE IV. 

- Continued. 

SAMPLE OUTPUT 




(a) 


' MN ARRAY <RhU*W- 
25.507595 

sub-mi 

26. 1C7347 

26.522294 

26.963462 

27.174644 

27.2C1535 

27.162536 

26.956273 


26.821924 

27.C283C3 

27.457889 

27.935325 

28.133250 

28. 1615CE 

28. 120715 

27.907177 

11 1 

26.538462 

27.582843 

2 8.347672 

28.937641 

29.10 5428 

25 .CEE 145 

25.069426 

28.951727 

25.099002 
26.28991 C 
„ 35.465356 

27.521081 

31.072922 

33.768741 

29.323885 
31.268500 
33. C45494 

30.027952 

31.024394 

31.902327 

30.074642 

30.616576 

30.882637 

25.525236 
3C. 137543 
25 . 5635C4 

29. 897213 
32.223495 
31.077170 

30.232534 


■ 1 

f MMU /FRAY <RHJ*h- 
16.971501 
! 53-571274 

55.5157C1 
54.498377 
47.28150C 

SUB-M ON UPPER 
37.7E7557 
54.158436 
56.160765 
54 .088461 
42.868861 

SURFACE) 

43. 857512 
54.195787 
56.063219 
53.399481 
7. 3233204 

47.694059 
54. 121918 
55.506120 
52.632174 

50.197552 

53.969568 

55.010693 

52.010118 

•1.845SC3 
52. £3*425 
54. £12426 
5 1 .663522 

52. 522299 
53.934666 
54. 840794 
52.421168 

53.588732 

54.563338 

54.742220 

52.255233 

11 i 

C UML ARRAY (RFO*W-SUB-M ON LOWER 
37.165581 33.425342 

| 15.943203 21.624109 

l 36.146715 37.31543C 

34.134425 34.211736 

L 35.15C435 34 . 16644 £ 

SURFACE) 

24.599400 
23.575663 
37.645579 
34.474140 
7. 5 8CC489 

20.346128 

25.699009 

37.362323 

34.625603 

18. 338b59 
27.890045 
36.700486 
34.322540 

17.635555 
3G.103C14 
3 5 .C 1 £C 16 
33.55EC13 

17.82C372 
32.292383 
34. 93626 £ 
32.302812 

18.639874 

34.372876 

34.317511 

34.168130 


11 



A PRAY ( RHU*W- SUB- ThE 


-35.624538 -35. 
-35.059493 -34. 
-30.427464 -29. 
-27.355055 -23. 
-12 .940090 -16. 


TA ) 


9C5356 

-4C.041 C67 

775723 

-34.638582 

084725 

-2 £• 989547 

802847 

-23.108671 

356022 

-18. 165348 


-39.991943 

-34.688226 

-29.231843 

-23.714139 

-19.045188 


-39.922701 

-34.758195 

-29.475374 

-24.209203 

-19.242202 


-35 

6EC27C 

-39 

-34 

E01C71 

-34 

-25 

6C2664 

-29 

-24 

3952C6 

-24 

-18 

.521674 

-18 


848412 -39.774967 
833264 -34. 90740 L 
676812 -29.923152 
407707 -24.706454 
543125 


11 


MXU X L PXL 

(H C00RU. VS. R KD* W- SUB- THE 1A ON UPPER SURFACE) 


-0 

-25.724336 

0. 442 3 5 76 E— 01 

0 .545968*6-01 

-3E. 356696 

0.5651698fc-Ul 

0 .6096762 £—0 1 

— 61 .064705 

0. 621 84186-01 

0.652922CE-01 

- 71.5C9989 

0.6617128E-01 

w 0 .684266 6 E— C 1 

- 113. 30783 

0.6 84 40 COE-01 


WXU 

MXU 

WXU 

MXU 

WXU 

-12.260539 
—44 .427 730 
-66.415494 
-73.904902 
-82.64 8045 

0 .4918970 E-C 1 
0 .58172576— C 1 
0.6330672E-C1 
0 .6698656E-C1 

-23-14*364 
-5C. 225415 
-6E. 043138 
-75.55C646 

0 . 52 2 6 6 0 76 -01 
0.5963962E-01 
0.64340326-01 
0.6774543E-01 

-31.450157 

-55.820843 

-69.741194 

-77.330525 


11 


(M 


MXL 

COURU. VS 
0 .4229608 
0 .55843*5 
0.61563C7 
0.6528332 


MXL 

. R FO* W- SUB— THE TA 
E-C 1 -fc. 6572541 
E— 0 1 -27.8U6565 

E— 0 1 — 45.4U4C54 

E— 0 1 -65.0E7313 


ON 


MXL 

LDwfcR SURFACE) 
0. 4 61 3 633E— 0 1 
0. 57586B6E-01 
0. 62614 71E-01 
C. 66053946-01 


MXL 

-13.911637 

-31.960470 

-50.232773 

-72.164698 


n*L 


0 .5133148E-0 1 
0 . 590 8518E-01 
0.6357468E-C1 
0.6677985E-C 1 


-15.51 63 CC 
-36.318843 

— 55 . 1 15CC5 

— 85.6354)6 


0.5379 02 IE-01 -23.899858 

0.60395 99E -01 -42.909056 

0.6445978E— 01 -60.541574 


11 



OF RHU*W 
.342670 
.142764 
.638791 
.12 7857 
.30 196 £ 


AT 


INIERIOR POINTS 
47.686801 
4 4. 04 4 Co9 
4C. 064095 
36.386610 
35.114754 


46.028316 
44.201437 
4C. 546362 
37.334982 
36. 162120 


48.243784 

44.538247 

41.132562 

38.262753 

36.403739 


48.293720 

44.717020 

41.423708 

38.607895 

36.161265 


46.274C24 
44.768126 
4 1 .5C24C? 
3£ .6 1 14 1C 
35.5E5 126 


48.225504 
44.767521 
41. 542083 
38.595071 
37.177965 


48.048815 

44.691640 

41.636492 

39.043757 


a See p. 20. 


33 


I 


II 


TABLE IV. - Continued. SAMPLE OUTPUT 


Descrip- 

tion 

(a) 


VELOCITIES AT INTERIOR HE SF PCINTS 


IA- 

1 

VELOCITY 
1J0.SC 
133. 2 C 

ANGLE COEG) 
-56.82 
-55. 7C 

VELOCITY 
131. 5C 
133.06 

ANGLE (DEG) 
-56.81 
-55.72 

VELOCITY 

132.49 

132.55 

ANGLE! DEG) 
-56.48 
-55.07 

VELOCITY 

133.11 

ANGLE (UEG ) 
-55.99 

VELOCI TY 
133.26 

ANGLE (DEG) 
-55.76 

IA- 

2 

VELOCITY 
122.26 
124. Ci 

ANGLE (DEG ) 
-52.58 
-51.02 

VELOC1 TY 
121. 98 
124.07 

ANGLE (DEG) 
-52.15 
-51.09 

VELOCITY 

122.43 

123.85 

ANGLE! DEG) 
-51.60 
-51.36 

VELOCITY 

123.41 

ANGLE (OEG ) 
-51.15 

VELOCI TY 
1 23. 93 

ANGLE (DEG) 
-51.01 

I A“ 

3 

VELOCITY 

113.03 

115.56 

ANGLE (DEG) 
-48.48 
-45. 5C 

VELOCI TY 
111.47 
H5.6e 

ANGLE (DEG) 
-46.52 
-45.59 

VELOCITY 

112.80 

115.95 

ANGLE! DEG) 
-45.64 
-45.95 

VELOCITY 

114.49 

ANGLE (DEG) 
-45. 2S 

VELOCI TY 
115.33 

ANGLE (DEG) 
-45.36 


ER AT ION NO. 

12 

MAXIMUM 

RELATIVE CHANGE IN DENSITY = 

0.4274E- 

03 






* 


SURFACE VELCCI TIES 

BASEC UN 

AXIAL 

COMPUNENTS 





* 


UPPER SURFACE 


* 


LOWER SURFACE 


H 

* 

VELOCITY 

ANGLE (OEG) SLRF. LENGTH 

RHC*W 

• 

VELOCITY 

ANGLE (OEG ) 

SURF. LENGTH 

RHG*W 

0 

* 

0 

SO.OC L 

16.972 

* 

0 

—SC .CO 

0 

37.170 

0. 20136-05 

* 

1C7.29 

1.8S 0.2126 E— 02 

37.808 

* 

94 • 5 6 I 

- 1.96 

C. 2126E-02 

33.454 

Q.4026E-G2 

* 

126.19 

1.96 C.4140E-02 

43.884 

* 

69.657 

-2.C4 

0.4140E-02 

24.614 

0.6039E-C2 

♦ 

138.02 

1.96 0.6154E-02 

47.722 

* 

57.869 

-2.C3 

C.6154E-C2 

20.358 

0.8052E-C2 

* 

147.59 

1.8S C. 8168E-02 

50.225 

* 

52 .445 

-1.95 

0.8168E-02 

18. 349 


surface velocities based un tangential cghfcnents 

UPPER SURFACE 

H VELOCITY ANGLE (DEG I RHC*W 


*24* VALUE CF RHO*W IS TOO LARGE* 


CECK NAME 
. ERROR 
CENSTY 
TASVEL 
2 CCP 


IFN OF CALL 
C0C25 
00G16 
00 524 

otoie 


ABS. LCC . 
1661 5 
IC132 
26113 
C3056 


-0 

0.4424E-C 1 
0.4919E-C 1 
0.5227E-C 1 
0 .5460E-C 1 


71.836 
17S.CS 
ISC. 04 
20C.E 7 
212.49 


so.cc 
-12.61 
-23.06 
— 30. 50 
-36.16 


25. 724 
56.165 
5S.095 
61.975 
65. G06 


14 


M 

0.4230E-C 1 
0 .4814E-C 1 
0.5133E— C 1 
0.5379E-C 1 
0 .5584E-C1 


LOWER SURFACE 


VELOCITY 
135.46 
114. C€ 
123.33 
131.96 
139.23 


ANGLE (DEG) 
-8. 87 
-22.23 
-29.54 
-34. 52 
-38. 66 


PHC*h 

43.442 

36.779 

39.579 

42.172 

44.316 


a See p. 20. 
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TABLE IV. - Concluded. SAMPLE OUTPUT 


Descrip- 
tion 
(a) _ 


BLADE SURFACE VELOCITIES 


bCm 

1 

I 


1 + 

1 + 
1 

C.C2C1-- 

1 

1 

1 

1 

I 


14 < 


LOO. 

♦ l* 

L 

1 

1 

1 

1 


150. 


35C . 
— 1 — 

1 


450. 

L 

1 

1 

l 


500. 

1 

1 

L 

1 

1 

1 

1 

1 

1 


550. 

1 

1 


C . 0401 

1- 

-+ 

-L 

-* 1 

i 

1 

1 

1 

1 

1 — 

— 

— 1 — 
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Error Conditions 


The error message is given first for each error condition. 

(1) MX, NBBI, NUSP, NLSP, NRSP, OR NINT IS TOO LARGE. If this message is 
printed, reduce the appropriate input values to the stated maximum value. 

(2) WTFL IS TOO LARGE AT UPSTREAM BOUNDARY is printed out if WTFL (w) is 
greater than the choking mass flow for the upstream boundary AH. If the "continue” 
control card is used (see p. 86), WTFL will be cut in half, and calculations will proceed. 
This allows the calculation of further useful information. 

(3) INU AND INL MUST BE LESS THAN 100 is printed if there are more than 100 in- 
tersections of horizontal mesh lines with either the upper or the lower blade surface. In 
this case NBBI should be decreased. 

(4) INU NOT EQUAL TO JU OR INL NOT EQUAL TO JL. The number of intersec- 
tions of horizontal mesh lines with either the upper or lower blade surface are counted by 
both COEF and TASVEL. If these counts do not agree, the above error message is 
printed out. This error is probably due to an error in the input. 

(5) THE NUMBER OF UNKNOWN MESH POINTS EXCEEDS 2500, A COARSER MESH 
MUST BE USED is printed if there are more than 2500 interior mesh points. The actual 
number of interior mesh points is given. Either MX or NBBI must be reduced. 

(6) VALUE OF RHO* W IS TOO LARGE is printed if the value of pW at some point 
is so large that there is no solution for the value of p and W. Decreasing WTFL (w) 
sufficiently eliminates this condition. However, it may be desired to continue the calcu- 
lations. If so, the "continue" control card (see p. 86) is used. This may permit an ap- 
proximate solution to be obtained, which would be valid at other points. In some cases 
the value of pW is reduced at the point in question during later iterations, resulting in a 
valid final solution. 

(7) OUT OF RANGE Z = x.xxx is printed if SPLINT is used for extrapolation. Also, 
the input and output for SPLINT is printed. SPLINT is normally used for interpolation, 
but may be used for extrapolation in some cases. Calculations proceed normally after 
this print out. 

(8) CAUTION -HB* RM(MXBO) LESS THAN RO LESS THAN HA MAY NOT GIVE 
CORRECT RESULTS is printed if the internally calculated values of HA and HB are such 
that 


HB* RM(MXBO) < RO< HA 

Decreasing NBBI or increasing the difference MXBO - MXBI alleviates this condition. 
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PROGRAM PROCEDURE 


The program is segmented into five main parts - the subroutines INPUT, COEF, 

SOR, SLAXVL, and TASVEL called by the main program 2DCP. In addition there are 
several other subroutines. All the subroutines and their relation are depicted in fig- 
ure 13. All information which must be transmitted between the five main subroutines is 
placed in COMMON. The program can handle up to 2500 mesh points on the IBM 2-7094- 
7044 direct coupled system with a 32 768 word core. 


2DCP l 


| TIME1 


INPUT^ | COEF 


| SPLN22| 


DEDERR 


| INTPLlj — 
\ BLDCR1 | 


ARERR 


Qor] 


SPLINT 


DEN STY \ 

_J 


SLAXVL | | TASVEL 1 


| VELOC 


I 


SPLINE 


I PLOTMY 
PISTUG 

~TZ 

| KHAR 


SORTXY 


BLDDE1 


Figure 13. - Logical relation of subroutines. 



Figure 14. - Arrangement for overlay, showing octal storage requirements. 
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To be able to handle 2500 mesh points an overlay arrangement is used as indicated in 
figure 14. All subroutines not shown are in the main link. The total program storage re- 
quirements is 74303^ of which 47106^ is in COMMON. The system storage require- 
ment for our computer is 2764^ and unused storage is 511^. If there is a storage prob- 
lem on the user’s computer, the maximum number of mesh points should be reduced. 

The first segment of the program is INPUT. This subroutine reads all input data 
cards, calculates basic constants and useful information, and calculates the blade coor- 
dinates on mesh lines. INPUT is called only once for each case. The next subroutine is 
COEF which calculates the entries of the matrix A and the vector k of equation (A7). 

These coefficients must be recalculated for each outer iteration. The subroutine SOR es- 
timates an optimum overrelaxation parameter £1 on the first call if it is not given as in- 
put. The same value of is used for each outer iteration. SOR then finds the linear 
solution to equation (A7) with fixed coefficients by successive overrelaxation. Then 
subroutine SLAXVL calculates the streamline locations and pW m and plots the stream- 
line locations if desired. Finally, the subroutine TASVEL calculates pW 0 , velocity 
magnitudes and direction, densities for next outer iteration, the surface velocities 
based on axial velocities and on tangential velocities, and plots the surface veloci- 
ties. 


Conventions Used in Program 

For convenience, a number of conventions are used in naming variables and assigning 
subscripts. First, several pairs of variables are spelled the same except for one letter, 
which is U in one case and L in the other. The U signifies the upper surface BC, and L 
the lower surface CF. Another practice is to use the letters I and O in a similar manner, 
where I refers to the inlet or the region ABGH, and O refers to the outlet or region 
CDEF. Thus, ALUO refers to the angle on the upper blade surface near the outlet, or 
near point C (fig. 10, p. 15). 

The variable I is used to number all the mesh points starting with I = 1 at A and pro- 
ceeding along the vertical mesh lines and moving to the right to the next line after the end 
of each vertical line and ending with I = NXN at the last mesh point near E. The mesh 
spacing in the m direction is labeled HA, and the spacing in the 9 direction is HB. 

The techniques used in the program and correspondence to the mathematical equa- 
tions are described briefly. Each subroutine is described separately, first the five seg- 
ments of the main program, followed by descriptions of each of the remaining subroutines. 
The various segments of the subroutines are labeled by comment cards, which generally 
correspond to the headings in the following descriptions. 
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Subroutine INPUT 


Input. - The first step is to read all input cards for a particular case. A detailed 
description of the input required is given in the section Instructions for Preparing Input 
(p. 14). All input data are given as the first output. 

Cal culat ion of cons tants and initialization . - After all input has been read in, the 
various constants needed in the program are calculated and certain quantities are ini- 
tialized. The input 0-coordinates for the lower blade surface (XSPL) are increased by 
PITCH to define the blade passage. Certain convergence tolerances are specified at this 
point. The arrays RM and BE of the quantities r and b at each vertical mesh line are 
calculated by subroutine SPLINT (using cubic spline interpolation) . Prerotation X is 
calculated by an interative procedure. Also, the useful quantities of upstream and down- 
stream free-stream velocity, maximum values of the mass flow parameter pW, relative 
critical velocity, and free-stream flow angle corrected to leading and trailing edges are 
all calculated and printed. All density arrays are initialized to p! n (RHOIP) . 

Calculation of mesh co o rdinates along boundary . - The 0-coordinates of boundaries 
BC and GF at each vertical grid line are calculated by BLDCR1 and stored in the arrays 
XU and XL. BLDCR1 requires as input the first and second derivatives at each spline 
point of the cubic spline curves describing the blade surface. These values are calcu- 
lated by SPLN22. The first and last points of the spline curves are determined by the 
angle of tangency to the leading- and trailing-edge radii. Therefore, these points are not 
specified as input, but are computed by this section of the program. 


Subroutine COEF 

The coefficients of u in equations (A 2) to (A 6) (elements of matrix A in eq. (A 7)) are 
computed at the same time as the constants (components of k in eq. (A7)). Between the 
blades it is necessary to compute values for hg and h^ at some mesh points adjacent to 
the boundary. These values are calculated by INTPL. Also it is necessary to calculate 
values for bg and b^ along the boundary . These values are calculated by SPLINT and 
stored in BEU and BEL. These arrays are ordered by increasing m-coordinate. This 
is not the same order as the order of mesh points (i.e. , increasing I subscript) which 
necessitates some juggling of the subscripts of the BEU and BEL arrays. The subscripts 
INL and INU increase with I, and ITP is the correct subscript for the BEU or BEL array. 
A similar situation holds with the RHOLT and RHOUT arrays for the values of Pg and 
p^ from the previous iteration. 

Near the trailing edge, a special situation may arise, as illustrated in figure 15. 

Here it should be noted that the blade intersects the mesh line twice between two adjacent 
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mesh points due to the small trailing-edge radius. This situation would not be detected 
by the program in the normal procedure, which leads to a large error in the velocity cal- 
culation at the next to the last vertical grid line on the lower blade surface. Therefore, 
a special check is made for the two vertical mesh lines involved at statements 120 and 
130, and if this situation occurs, the proper values of H3 or H4 are calculated. Also, 
the number of the horizontal mesh line is stored in IBTE2 and IBTE1, and this informa- 
tion is used in calculating the tangential velocity components. If the values of HA and 
HB are such that two horizontal mesh lines could intersect the trailing-edge radius, the 
message for error condition (8) (p. 36) is printed out. 


Subroutine SOR 


Estimation of value of optimum over relaxation factor. - If a value of W> 1 is given 
as input, it is used for the overrelaxation factor. Otherwise a value is estimated in the 
first iteration by using equation (B3) of reference 5 to estimate the value of p 2 (B) =p(L 1 ) 
and equation (Bl) of reference 5 to obtain the corresponding value of W (see appendix B 
of ref. 5). Equation (A 8) is used to calculate u m+1 from u m for equation (B 3), with 
oj = 1 and k = 0. To start, u° = 1 for all i. Equation (A8) becomes 


i-1 

n 


-r 1 - - Z 

vr 1 - £ vr 

(7) 

i-i 

j=i+l 



In the program, i is replaced by I directly. For each i, there are only four values of j 
for which a.- is nonzero, which are the negative values of the coefficients A(1, 1), 

A (I, 2), A(I, 3), and A(I, 4). The value of j is determined by the index of the proper 
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neighboring point. These indexes are named II, 12, 13, and 14, and are defined so u™ ^ 

has the coefficient A(1, 1), and similarly for the other coefficients. After the values of 

the indexes are computed, equation (7) is used to compute u m+ * from u m . Then, the 

minimum and maximum values of the ratio u [ n+ / uP are calculated and given the names 

LMIN and LMAX, respectively. After convergence, the optimum value of the overrelaxa- 

2 

tion factor Cl can be calculated from equation (Bl) of reference 5, since p (B) = LMAX. 

Calculation of initial solution estimate . - For the first outer iteration an initial solu- 
tion estimate must be made. This is done by assuming u = 0 along ABCD and u = 1 
along EFGH and assuming linear variation along vertical mesh lines. On subsequent iter- 
ations, the previous solution is used as the initial solution estimate. 

Solution of matrix e quation by SOR . - With a value of Cl either as input or estimated 
by the program, equation (A8) can be used iteratively to calculate a sequence ju m j that 
will converge rapidly to a solution of equation (A7) . The indexes i and j and the cor- 
respondence of a.j and uP to the program variables is the same as described pre- 
viously for estimating the optimum over relaxation factor. During each iteration the max- 
imum change of the stream function is calculated. When this maximum change is reduced 
below TOLER, set equal to 10 , the iteration is stopped, and the current estimate of the 

stream function is accepted as the solution. 


Subroutine SLAXVL 

Calculation of streamline locations and 9u /90 . - Along most vertical mesh lines, the 
stream function is a one-to-one function of the distance in the 0-direction. Therefore, 
the 0-coordinate is considered to be a function of the value of the stream function, and 
the value of 0 at a given value of the stream function can be obtained by cubic spline in- 
terpolation (SPLINT). At the same time, 9u/9 0 is computed along the same mesh line, 
estimating the derivative at each mesh point by use of the cubic spline (SPLINE). The 
derivative 9u/90 at unknown mesh points is stored in the array WM and 9u/90 along the 
blade surfaces is stored in WMU and WML. These calculations are performed in three 
sections, as noted by comment cards: (1) upstream, (2) between the blades, and 
(3) downstream. 

Plotting of streamlines . - The streamlines can be plotted to give a rough idea of 
their locations. These locations are particularly helpful in disclosing quickly any errors 
of input. The plotting printout is done by PLOT MY, which, with the necessary further 
subroutines PISTUG and KHAR, is described completely with FORTRAN IV listing in ref- 
erence 12. The plotting can be omitted by removing statements following statement 420 
up to and including statement 470. 

Calculation of pW m . - The product pW m is calculated by multiplying 9u/90 by 
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w/br (using eq. (3)). The values of pW m are stored in WM, WMU, and WML. 


Subroutine TASVEL 

Calculation of 3u/3m. - The tangential velocity component is calculated from 
(3u/3m) = (-bpW^/w) by considering each horizontal mesh line. The fact that the various 
horizontal mesh lines start and end at various places complicates this process. To sim- 
plify the procedure, upstream and downstream ends of each mesh line are considered 
separately. At the upstream end, there are three possibilities: (1) the line starts at 
AH(IA = 1) , (2) the line starts on lower surface of blade, or (3) the line starts on the up- 
per surface of the blade. Similarly, at the downstream end there are three possibilities. 
For convenience, case numbers are assigned to the various possibilities as follows: 


Case 

Starts 

Ends 

1 

AH 

Upper surface 

2 

AH 

DE 

3 

AH 

Lower surface 

4 

Lower surface 

DE or lower surface 

5 

Upper surface 

DE or lower surface 


As mentioned in the description of COEF, a special situation often arises where a 
horizontal mesh line is intersected twice between two adjacent mesh points, as illustrated 
in figure 15 (p. 40). When this occurs, the index of the horizontal mesh line is stored in 
IBTE2 and IBTE1. Under cases 2, 4, or 5, if IB = IBTE2, the special case arises, and 
previously calculated information is used for this line to the left of the trailing edge. 

After all other tangential velocities have been calculated, tangential velocities are calcu- 
lated for the remainder of this line (IB = IBTE1) to the right of the trailing edge. 

TASVEL calculates the necessary information about the two end points of each hori- 
zontal mesh line by using INTPL to calculate the mesh spacing at the end points. Then 
VELOC calculates 3u/3m at each point along the line by using SPLINE to calculate the 
actual derivatives. The derivatives, 3u/3m at unknown mesh points are stored in the 
array WX, and 3u/3m along the blade surfaces is stored in WXU and WXL, with the 
corresponding m -coordinates stored in MXU and MXL. The values of WXL and MXL 
are rearranged in increasing order of MXL by SORTXY. 

Calculation of pW g . - The product pW 0 is calculated by multiplying 3u/3m by 

-w/b (using eq. (2)). The values of pW 0 are stored in WX, WXU, and WXL for interior 
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mesh points, upper surface of blade, and lower surface of blade, respectively. 

Calculation of mass flow parameter pW and angles at interior points . - At each in- 
terior point, pW is calculated by pW = ^(p w m ^ 2 + ^pW^ 2 , and the angle 0 is calcu- 
lated by tan 0 = pW 0 /pW m . 

Calculation of velocity and density at each in terior point . - A value of pW deter- 
mines a unique subsonic velocity W and corresponding density p. These are calculated 
at each interior point by subroutine DENSTY. Also the relative change in p at each 
point from the previous iteration is calculated, and if the maximum relative change in 
density is less than 0.001, the outer iteration is considered to have converged sufficiently 
and the calculation terminates after the final printouts. 

The derivative 9u/9m at unknown mesh points is stored in the array WX. 

The derivative 3u/9m along the blade surfaces is stored in WXU and WXL, with the 
corresponding m-coordinates stored in MXU and MXL. 

Calculation of surface velocity based on pW m< - At each vertical mesh line pW at 

the blade surface is calculated by 


pW = 


P W m 
cos 0 



The derivative de/dm of the blade surface at each vertical mesh line is computed by 
BLDCR1 at the same time the blade coordinates XU and XL are computed, and is stored 
in DXDZU and DXDZL. The surface velocity is then calculated from the value of pW by 
DENSTY. The surface velocity based on pW m is more accurate at small values of 0 

and would not be expected to be accurate for 1 0 J > 60°. The blade surface length is cal- 
culated for convenience using equation (B17). 

Calculation of surface velocity based on pW 0 . - At each horizontal mesh line pW at 

the blade surface is calculated by 



The derivative de/dm of the blade surface at each horizontal mesh line is computed by 
BLDDE1 and is stored in DTDMU and DTDML. The surface velocity is then calculated 
from the value of pW by DENSTY. The surface velocity based on pW 0 is more ac- 
curate when 1 0 1 is close to 90° and would not be expected to be accurate for 1 0 1 < 30°. 
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Plotting of surface velocities. - If desired, the surface velocities are plotted using a 
printer plotter. The velocities are plotted using different symbols for upper and lower 
surface and for velocities based on meridional components or on tangential components. 
Velocities based on meridional velocity components are plotted if | /3 1 < 60° and velocities 
based on tangential velocity components are plotted if j /3 1 > 30°. Plotting is done by 
PLOTMY, which is described in reference 12. 


Internal Variables for INPUT, COEF, SOR, SLAXVL, and TASVEL 


A 

array of coefficients of u which are elements of matrix A 
in eq. (A7) 

A12, A34 

a 12’ a 34 in eq ‘ 

AA 

temporary storage 

AAA 

array used for temporary storage 

An 

a Q in eq. (A2) 

B 

temporary storage 

B12, B34 

b 12> b 34 in eq. (A 2) 

BE3, BE4 

b 3 , b 4 in eq. (A2) 

BE 

array of values of b at vertical mesh lines 

BEL(BEU) 

array of values of b at horizontal mesh lines on lower (upper) 
blade surface 

BETA 

array of values of /3 at interior mesh points 

BETA L(BETAU) 

array of values of /3 on lower (upper) blade surface 

BTAIN 

free-stream angle j 3 at blade leading edge based on /3. n , 
calculated by eq. (B13) 

BTAOUT 

free-stream angle /3 at blade trailing edge based on @ ou {., 
calculated by eq. (B15) 

CASE 

number (integer) of case in calculating tangential velocity 
components 

CHANGE 

change in value of stream function at a particular point when using 
SOR iteration 

CP 

C P 

CPTIP 

2 c p T in 
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DELINT 


DTDML(DTDMU) 

DTLR 


DX 

DXDZL(DXDZU) 

DZ 

EML(EMU) 

ERROR 

EXPON 

FIRST 

HI, H2, H3, H4 
HA 

HAMRO 

HB 

HA1 

HA1TE 

HAN 

HANTE 

HL 

HU 


increment of stream function for which streamline locations are to 
be calculated 

array of d0/dm at horizontal mesh lines for lower (upper) blade 
surface 

tolerance for mesh points near boundary mesh point (If a mesh 
point is closer than DTLR to the boundary, the boundary is con- 
sidered to go through the mesh point. The program uses 
DTLR = 0.001 HB.) 

0 plotting increment for streamline plot 

array of values of slope of lower (upper) blade surface at each 
vertical mesh line 

m plotting increment for streamline plot 

array of second derivatives of spline curve at each spline point 
for lower (upper) blade surface, calculated by SPLN22 

maximum absolute value of the change in u for an overrelaxation 
iteration 

i/(y - 1) 

value (integer) of I at lowest mesh point for given vertical mesh 
line 

hj, h 2 , h 3 , h 4 (see fig. 18 in appendix A) 

basic mesh space in meridional (m) direction 

HA-RO 

basic mesh space in blade -to -blade (0) direction 

length of first mesh space along horizontal mesh line 

HA1 for special case shown in fig. 15 for line segment to right of 
trailing edge 

length of last mesh space along horizontal mesh line 

HAN for special case shown in fig. 15 for line segment to left of 
trailing edge 

0 distance between EF on boundary and first mesh line below 
(fig. 12) 

0 distance between CD on boundary and first mesh line above 
(fig. 12) 
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I 


11(12,13,14) 

IA 

IB 

IBDL 

IBDU 

IBTE1,IBTE2 

IA1 

IAN 

IL(IU) 

INL(INU) 

ITE2 

ITER 

ITERA 

ITP 

J, JI, JB 
JL(JU) 

JUM1 

K 

K1,K2,K3,K4,K5 

KK1 


KN 

KKK 

LAMBDA 


index of mesh point 

index of mesh point located at 1 (2,3,4) in fig. 18 with I at 0 

index of mesh line in meridional (m) direction 

index of mesh line in blade -to -blade ( 9 ) direction 

difference between NL for a vertical mesh line and the next one 

difference between NU for a vertical mesh line and the previous 
one 

indexes of special mesh lines shown in fig. 15 

index of first mesh point along horizontal mesh line 

index of last mesh point along horizontal mesh line 

array of indexes of highest (lowest) mesh point for each vertical 
mesh line 

index counting number of intersections of horizontal mesh lines 
with lower (upper) surface 

index of mesh point indicated in fig. 15 

outer iteration number 

temporary storage of outer iteration number 
temporary index 
temporary indexes 

number of points where horizontal mesh line intersects lower 
(upper) blade surface 

JU-1 

array (real) of constants that is vector k in eq. (A7) 

code variables (real) used in determining values of coefficients 
A (I, J) and constants K(I) 

code to specify whether first point of horizontal mesh line is on 
AH(KK1 = 0) or upper blade surface (KK1 = 0) or lower blade 
surface (KK1 = 1) 

same as KK1, but for last point 

array containing information used in plotting subroutine PLOTMY 
X 
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LAST 


LMAX 
LMIN 
MPL 
MXBIM1 
MXBIP1 
MX BO Ml 
MX BO PI 
MXL(MXU) 

NBB 

NBBO 

NBUO 

NCH 

NL 

NP1(NP2) 

NP3(NP4) 

NSP 

NU 


NUTEMP 

P 

PITCH 

RATIO 


value of I at highest mesh point for given vertical mesh line 
upper bound (real) for p(L^) from eq. (B2) of ref. 5 
lower bound (real) for p(L^) from eq. (B2) of ref. 5 
array of m -coordinates of vertical mesh lines 
MXBI - 1 
MXBO + 1 
MX BO - 1 
MXBO + 1 

array of m -coordinates of intersections of horizontal mesh lines 
with lower (upper) blade surface 

number of mesh points along vertical mesh line 

number of mesh lines above mesh line AB for first mesh line 
below EF (may be negative) 

number of mesh lines above mesh line AB for line CD (usually 
negative, unless STGR is positive) 

number of vertical mesh lines in length of blade 

array of number of mesh points on vertical mesh line above line 
AB (may be negative) 

number of plotted upper (lower) blade surface velocities based on 
meridional components 

number of plotted upper (lower) blade surface velocities based on 
tangential velocities 

number of mesh points plus boundary points along vertical mesh 
line 

array; on vertical mesh line LA, the mesh point nearest the upper 
blade surface is NU(IA) mesh points above line AB (NU(IA) may 
be negative) 

temporary storage of NULAKI 

array of input information for plotting subroutine PLOTMY 
2 tt/NBL 

value of ur i+ ^/u| n for use in eqs. (B2) and (B3) of ref. 5 
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RELER 

maximum relative change in density at interior mesh points, be- 
tween two outer iterations 

RHOl, RH02, 
RH03, RH04 

pi, p2, p3, p4 in eq. (A2) 

RHO 

array of densities p at interior mesh points 

RHOL(RHOU) 

array of densities p at vertical mesh lines for lower (upper) 
blade surface 

RHOLT(RHOUT) 

array of densities p at hozizontal mesh lines for lower (upper) 
blade surface 

RHONEW 

newly calculated estimate of p 

RHOT 

temporary storage of a value of p 

RHOVI 

(O^in 

RHOWMI 

maximum value of pW along AH 

RHOWMO 

maximum value pW along DE 

RM 

array of values of r at each vertical mesh line 

RML(RMU) 

array of values of r at the intersections of horizontal mesh lines 
with the lower (upper) blade surface 

SAL 

array of values of sin a - (dr/dm) at each vertical mesh line 

SL 

array of streamline coordinates for input data to the plotting sub- 
routine PLOTMY 

SLLPE(SLUPE) 

array of slopes of spline curve at each spline point for lower 
(upper) blade surface, calculated by SPLN22 

SRW 

code (integer) variable that will cause certain subroutines to 
write out data useful for debugging: 

SRW = 13 SPLINE will write input and output data 

= 16 SPLINT will write input and output data 

= 18 SPLN22 will write input and output data 

= 19 BLDCRlwill write out blade coordinates and 

slopes at each mesh line 

= 19 INTPLlwill write out pertinent data for each 

iteration 

= 20 BLDDElwill write out m-coordinates and slopes 

TANTH 

tan /3 at unknown mesh points 

TANTHL(TANTHU) 

tan (3 along lower (upper) blade surface 
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TBI 

tan ^in 

TBO 

tan ^out 

TGROG 

2yR/(y + 1) 

TOLER 

when the maximum absolute value of the change in the stream 

function is less than TOLER, the SOR iteration is considered 

0 

converged; value used in the program is 10" 

TPP 

rpf t 

TTIP 

T / T in 

TWL 

2u)X 

TWLMR 

2 - (wr) 2 

TWW 

2co/w 

U 

array of estimated values of stream function or of eigenvector 
associated with spectral radius of Lj, p(Lj) as estimated by 
power method (ref. 5) 

UINT 

array of values of stream function for which it is desired to ob- 
tain interpolated values of 9 -coordinate 

UNEW 

new value of eigenvector estimate at single point, as calculated 
eq. (7) 

UNI 

(3u/ar?) in (eqs. (5) and (A3)) 

UNO 

(0u/0rj) out (eqs. (6) and (A 4)) 

USP 

array of values of stream function along vertical mesh line, 
including boundary points 

V 

array of relative velocities W at unknown mesh points, also 
used for storing values of pW 

VEL 

temporary storage, W 

VI 

W in 

VO 

W out 

VTOL 

in calculating W from pW by eq. (B6), the procedure is con- 
sidered converged when the relative change in W is less than 

_4 

VTOL; program uses a value of 10 

WCRI 

W£ r at B, fig. 4 

WCRO 

W” r atC, fig. 4 


49 



WL(WU) 

WM 

WMAX 

WMAX1 

WMIN 

WML(WMU) 

WR 

WX 

WXL(WXU) 

XI 

XBB 

XDOWN 

XFACT 

XINT 

XL(XU) 

XMAX 

XMIN 

YACROS 

ZINT 

ZFACT 

ZMIN 


array of velocities along lower (upper) blade surface 
array of pW m at interior mesh points 

upper bound for optimum from eqs. (Bl) and (B2) of ref. 5 
temporary storage for WMAX 

lower bound for optimum £1 from eqs. (Bl) and (B2) of ref. 5 

array of pW m where vertical mesh lines intersect lower (upper) 
blade surface 

tolerance specified for the calculation of overrelaxation factor n, 

-5 

program uses 10 

array of pWg at interior mesh points 

array of pW 0 where horizontal mesh lines intersect lower (upper) 
blade surface 

value of 0 for which INTPL1 is to compute H3 or H4 
array of 0-coordinates associated with array USP 
array of m-coordinates where surface velocities are plotted 
scaling exponent for streamline plot 

array of interpolated 0 -coordinates calculated by SPLINT and 
corresponding to array UINT 

array of 0 -coordinates of lower (upper) surface of blade at each 
vertical mesh line 

maximum value of 0 in streamline plot 
minimum value of 0 in streamline plot 
array of surface velocities plotted 
argument of BLDCR1 , not used in main program 
scaling exponent for streamline plot 
minimum value of m for streamline plot 
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MAIN PROGRAM 

CUMMUN SR W * ITER 
SRW = C 

CALL TIMEl(Tl) 

10 CALL INPUT 
20 CALL CClEF(NXN) 

IF( NAN .GT .2500 ) GO TU 1C 

CALL SOR 

CALL TIMEKT2) 

TIME = { T2-T1J/36C0. 

WRITE (6, 100C ) TIME 
CALL SLAXVL 
CALL TASVEL 
CALL T1ME1IT2) 

TIME = (T2— Ti)/36GC. 

WRITE (6, 1000) TIME 
I FI ITER.EQ .0 ) GU TO 10 
GO TO 20 

1000 FORMAT (1FL.7HTIME = ,F7.A,5H MIN. > 
END 


SUBROUTINE INPUT 

IA IS AXIAL INDEX 

IB IS BL ADE-TO-BLAOE INDEX 

I IS OVERALL INDEX 

FA IS BASIC AXIAL INCREMENT 

FB IS BASIC BLADE- T0-8LADE INCREMENT 

TFE FOLLOWING CODE IS USED FOR OUTPUT OPTIONS 

1 - LAST ITERATION ONLY IAFTER CONVERGENCE) 

2 - FIRST AND LAST ITERATION 

3 - ALL ITERATIONS 

REAL K,K1,K2,K3,KA,K5,LMAX,LMI N, LA MBDA ,MU , ML , MR , MXU ,MXL,MPL 
INTEGER SRW,FIRST,CASE f BLDATA,ERPRT ,STRFN »SLC RD * SLPLT , ARPRT ,SU RV EL 
COMMON SRW, ITER , GAM , A R , T I P , RHOI P »WTFL , OMEGA , LAM EDA , CP , EX PON, P IT CH, 

1 CHORD, STGR.BETAI ,B E TA 0 ,D T LR ,RI ,ALUI , ALLI ,RC, ALUC , ALLO, 

2 MX B I.MXBO ,M X , NBB 1 , NUSP , NL SP , NR SP , M NT , VTCL , 

3 BLDATA ,NULAKI ,ERPRT ,STRF N ,SLCRD , ARPRT ,1 NTVEL ,SURVEL, 

A MU ( SO) ,XSPL(50) ,ML( 50) ,XSPL<50) ,MR(50) , R MS P 1 5 0 ) ,tJESP(50), 

5 W, WR, TOLER ,B0A ,BDD , U ( 2 500 ) ,A (2500 ,A) ,K(2500 ) , PH0I2500 ) , 

6 DXDZU( ICO ) ,DXDZL ( 100) ,SLUPE (50) ,EMU (50) ,SLLPE (50 ) ,EML(50 ) , 

7 RM( 100 ),BE( 100) , SA L ( 100) , XU ( 1 00 ) , XL ( 100) , RMU (100 ) , RML( 100 ) , 

8 NU( ICO ),NL( 1CC),UINT(11) ,XINT(11) ,MPL(100) , 

9 HA, FB,NXN,MXBIM1,MXB0P1 ,JU,JL,HU,HL,NBBC,NBUC , NCH , 

1 IBT El, IBTE2, I TE2.HA 1TE ,HANTE ,UNI , UNC , TWW , I T ERA , 

2 KHOUI ICG) ,R HOL l ICO) ,RHQOT (ICO) ,RHCLT (100) , BEU (100 ) , BEL ( 100 ) , 

3 AAA(IOC) 
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C THE FOLLCWING VARIABLES ARE ALL IN COMMON BY EQUIVALENCE 

DIMENSION WM{ 2500) ,WX(250C) , V(2500) , BE TA (2500 ) ,S L ( 1 100 ) 
EQUIVALENCE (Al j ,(A(1 ,2) ,MX(I) I , (All ,3) ,V(I) ) . 

1 ( A ( 1» 4 I » BETA! 1)1 , (K< 1401) ,SL(1) ) 

1C WRITE (6, 599 ) 

ITER = 0 

READ ( 5, 10 10 ) GAM.AR , TI P ,RHOl P , WTFL , OMEGA, W 
WRITE (6,1100) 

WRIT £(6,1C20) GAM.AR , TI P ,RHOI P , WTF L , CMEGA , W 
WRITE (6,1110) 

READ (5,1010) CHORD, STGR , BE TAI ,BETAC 
WRITE( 6, IC20) C HURD, STGR, BETAI .BETAC 
WRITE (6, 1120) 

READ (5,1010) RI.ALUI ,AL L I ,R0 , A LUO , AL LC 
W R I T E ( 6,1020 ) RI.ALUI , AL LI ,RC , A LUO , ALLC 
WRITE (6, 1130) 

READ (5,1000) MXBI,MXBO,MX,NBBI ,NUSP,NLSP,NRSP,N8L,NINT 
WRITE! 6, 1C0C ) MXBI , M XBO , MX , NB8I , NUSP , NL SP , NRS P , NBL , M NT 
I F( MX .C-T. 100 .OR .NB8I .GT. 50. OR. NUSP. GT. 50. OR. NLSP.GT . 50 . CR. NRS P .GT . 
1 50.0R.NINT.GT.10) GO TO 3C0 

WRITE (6,1140 

READ (5,1010) ( MU(IA) ,IA=1,NLSP) 

WRITE! 6, 1020) ( M U ( I A ) , I A =1 ,NL SP) 

WRITE (6, 1150) 

READ (5,1010) (XSPU(IA) ,IA=1 ,NLSP) 

WRITE! 6, 1020) (XSPU(IA) ,IA=1 ,NLSP) 

WRITE (6,1160 

READ (5,1010) ( ML( IA ) ,IA=1 ,NLSP) 

WRITE! 6, 1020 ) ( ML ( I A ) , I A =1 ,NL SP) 

WRITE (6,1170) 

READ (5,1010 (XSPL(IA) ,IA=1 ,NL$P) 

WRITE! 6, 1020 > ( XSPL ( I A ) , I A =1 ,NLSP> 

WRITE (6, 118C ) 

READ (5,1010) ( MR ( IA) ,IA=1 ,NRSP) 

WRITE! 6, 1C20 ) ( MR ( IA) ,IA=1 ,NRSP) 

WRITE (6,1190 

READ (5,1010) (RMSP(IA) ,IA=1 ,NRSP) 

WRITE( 6, 1020 ) (RMSP(IA) , I A =1 ,NRSP) 

WRITE (6,1200 

READ (5,1010) (BESP(IA) , I A =1 »NRSP) 

WR 1TE( 6, 1C2G ) ( BE SP ( I A ) , I A =1 ,NR SP ) 

WRITE (6,1210 

READ (5, 1CCC ) BLDATA .NULAKI ,ERPRT, S TRF N , SLCRD , ARPRT , I NT V EL , 

1 SURVEL 

WRI T£( 6, 1CC5 ) BLOATA , NULAKI ,ERPR T , ST RF N , SLCRD , AR PRT , I NT VEL , 

1 SURVEL 

END UF INPUT, CALCULATION OF CONSTANTS AND INITIALIZATION 

PITCH = 2.*3.1415927/FLOAT(NBL) 

CO 15 IA = l.NLSP 
15 XSPL(IA) = XSPL ( IA ) +P I TCH 
FA = CFORD/FLQAT (M XBO— MXB I ) 

HB = P ITCH /FLOAT (NBBI ) 

DTLR = H6/10C0. 

A = ( P ITCF+STGR )/HB 
B = iTCR/FB 
NBBO = A- S IGN ( D TLR »A ) 

IF( A.L T.OTLR ) NBBO = NBBO-1 
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NflUO = 8-SIGN(DTLR,B ) 

IFIB.LT.— DTLR) NBGO = NBLO-1 

HU * FLO A T I NBUO + 1 )*HB— STGR 

hL = P ITC h*S T GR— FLQA T { NBBO ) *H8 

BETAI = BETA 1/57. 29 577 

BETAU = 8ETAO/57. 29577 

ALUl=ALUI/57. 29577 

ALLI=ALLI/57. 29577 

ALUU=ALUO/57. 29577 

ALLO=ALLO/57. 29577 

NCH = MXBO-MXBI+1 

MXBIM1 = MXBI-1 

MX8UP1 = MXBO+1 

IBTE1 = 1CCC 

IBTE2 = 1CCC 

WR = .COCCI 

TOLER = . CCGCC 1 

VTOL = .OCC1 

DO 20 IA= 1 »M X 

20 MPLIIA) = FLOAT! IA-MXBI )*HA 
CP = AR/( GAM- 1. )*GAM 
EXPON = 1 . / ( GAM— 1 • ) 

TWW = 2 .*OMEGA /WTFL 

CP T IP = 2.*CP*TIP 

TGROG = 2 .*GAM*AR /(GAM+l.J 

CALL SPL INTI MR ,RMSP ,NRSP ,MPL ,MX,RM, SAL) 

CALL SPLINT I MR , BE SP , NRSP , MPL , MX ,BE , A AA) 

TBI = SIN(BETAI)/ CO S ( BE TAI ) 

T BO = SINIBETAOJ/COSIBETAO) 

UNI = T8I/PITCH/RMI1) 

UNO = -TBO/P ITCH/RMI MX) 

ALU 1= S IN ( ALU I ) /CO S ( AL U1 ) /RM (1 ) 

ALL I = S IN ( ALL I ) /CO S ( A LL I ) /RM (1 ) 

AlUO = S INI ALUO ) /COS(ALLO) /RMIMX) 

ALLO = S INI ALLO ) /CO S I A LLO ) /RMIMX) 

C CALCULATE LAMBDA AND VI 
25 RHOT = R HO IP 

RHOVI = WTFL /BEI D/PITCH/CUSIBETAI ) /RMI) 

30 V I - RFOV I/RHOT 

LAMBDA = RMI 1 )* I V I* S I N (BE TA I ) + Cl MEG A * R M 1 ) ) 

T T IP = l.-I VI**2+2.*OMEGA*LAMBDA-(UMEGA*RM(l) )**2)/CPTIP 
IF! TT IP.LE.O.) GU TO 35 
RHONEW = RhO IP* T T IP**EXPON 

I FI AB S (RHONE W-R HOT ) /RHUI P.LT..C00UO1) GC TC AO 
RHOT = RhUNEW 
GO TO 20 

35 WTFL = WTFL/2. 

CALL ARERRI 39HWTFL IS TOO LARGE AT UPSTREAM BOUNDARY*) 
WRITE! t, 1A0G ) W1FL 
GO TO 25 

C CALCULATE MAXIMUM VALUE FOR RHO*W 
4C VI = RFOV I/RHONEH 

LAMBDA = RMI 1)*( VI*SIN(BETAI )+OMEGA*RM (1 ) ) 

TWL = 2.*OMEGA*LAMBDA 

AA = I TWL— IOMEGA*RM I 1 ) )**2 ) /CPTI P 

TPP = T IP * I 1 .— A A ) 

B = TGROG+TPP 
T T IP = I.-B/CPTIP-AA 
RHOT = RHOIP*TTIP**EXPON 


53 



c> r> o 


MIIHIIIII 


RHOWMI = RHO T* SUR T t B ) 

A A = { TWL-(OMEGA*RM(MX) )**2)/CPTIP 
TPP = TIP*( l.-AA ) 

8 = TCROG*TPP 
T T IP = l.-B/CPTIP-AA 
RHOT = RhO IP*TTIP**£XPON 
RhOWMO = RHO T* SQR TIB ) 

C CALCULATE VO AND ^-CRITICAL 

RhQVI = WTFL/BE(MX)/PITCH/COS(BETAO> /RM(MX) 

RHUT = RHO IP 

TWLMR = TRL- (OMEGA*RM ( MX ) ) **2 

CALL CENSTY (RHUVI ,RHOT, VO ,TwLMR ,C PTI P ,EXPCN, RHCIP »GAM, AR,T IP, 

1 VTOL) 

WCR1 = SURT( TGROG*TiP* (1.- ( TWL- (OMEGA*RM (MXBI ) ) **2 )/CPT IP) ) 
wCRO = SQRTI T GROG* TIP* (l.-(TWL— (OMEGA* RM(M XBO) )**2 )/CPT IP) ) 

C CALCULATE BETA CORRECTED TO BLADE L£ CR TE 

TBI = (LAMBDA-OMEGA*RM(MXBI )**2) ZW TF L*RHGNEW*BE ( MXB I ) *PITCH 
BTAIN = AT AN ( TBI ) * 5 7 .29577 

TBU = (TBU/BEIMX ) +OMEGA* (R M ( MX ) * *2 - RM ( MXBC) **2 ) *PHOT/wTFL* 

1 P ITCH )*BE(MXBO ) 

WRITE ^ ( 6, 1220 ) T LAMBOA VI »RHOwMI ,BTAI K ,I*CRI » VC , RHCWMC, BT AOUT , WCRO 

WRITE (6,1230 HA ,HB ,HU,HL ,PI TCH ,NBBC ,NBUU 

IF ( HB*RM (MX60 ) .LT .RO.AKD.RO. LT.HA ) uRITE (6,1240) 

DO 50 IA= 1 , ICC 
RHUU( IA ) = R HO IP 
R HOL ( I A ) = RHO IP 
RhuUK IA ) = RHO IP 
RHOL T ( IA ) = RhU IP 
BEU( IA ) = C . 

BEL ( IA ) = C. 

5C CONTINUE 

CO 60 1=1 , 25CC 
6C RHO( I ) = RHO IP 

IF( BL CATA .GT .C )WR ITE (6, 12 50 (PPL (I A) ,RM ( I A) , S AL U A ) , BE I I A ) , 

1 AAA( IA ), IA=1,MXJ 

CALCULATE MESH CuOROINAfES ON BOUNDARY 


AA = R M ( M XB I )*ALUI 

M U ( 1 ) = K I*( l.-AA/SQRT(l.+AA**2) ) 

XSPU(l) = R I/SQRT( l.+AA**2) /RM(MXBI ) 

AA = RMIMXBI >*ALLI 

ML ( 1 ) = R l*( l.+AA/SQKT(l.+AA**2) ) 

XSPL(l) = -R I/SUKK l.+AA**2) /RM(MXBI) + PI TCH 
AA = RM(M XBO )*ALUU 

MU(NUSP) = CHORD-RU* ( l.+AA/SGRT (1. + AA**2) ) 

XSPU(NUSP) - RU/SURT ( 1. +AA**2 ) /RM ( MXBC) +STGR 
AA = RM ( M XBO )*ALLU 

ML(NLSP) = C HURD-RG* ( 1 .—A A /SORT (i.+AA**2) ) 

XSPL(NLSP) = -RU/SURT( l.+AA**2) /KM(MXBC) +STGR+PITCH 
CALL SPLN22(MU»XSPU,ALUI , A L LO , NUSP , SL U FE , E MU) 

CALL SPUN 22(ML ,XSPL,ALLI ,ALLC, NLSP, SLLFE ,EML) 

CALL BLDCR1 ( MU , XSP L, SLUPE »E MU , NUSP »RI , ALUI ,R0 , A LU C , CHORD ,ST GR, 
1 PITCH, l.,X0,NCH,ZINT,M>3I ,D XDZU , RM ( MXB l ) ,RM(MX8CJ l) 

CALL BLDCR1 (ML , X SPL , SLL PE ,EML,NLSP,RI ,ALLI , RQ,ALLU, CHORD, ST GR, 
1 PITCh, — l.,XL ,NCH , Z I NT ,M XB I ,DXUZL,RM(MXBI) ,RM(MXBO) 

IF(BLCATA.LE.C) GO TO 65 
WRITE! t , 126C ) 
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WRITE! t, 1270 )(MU( IA),XSPU(IA) ,SLuPE (1A> ,EMU(IA) ,IA=i,NUSP) 
wR Iiet 6, 12 6 5 ) 

WRITE! <,12 70 ) (ML (IA),XSPL!IA) , SLLPE (1A) ,E ML { 1 A ) , I A= 1 , NLS P ) 
WRITE! 6, 1280 } <MPL! I A ) , XU < I A ) ,0 XD2U ( I A ) , XL ( 1 A> ,0X0ZL(I A) , 

1 IA = MXB I,MX6C ) 

6 5 CON TIN IE 
C CALCULATE NL,NL 

CO 70 1A = 1 * M XB 1M 1 
NL ( IA ) = N88 1- 1 
NO ( IA ) = C 
XU( IA ) = C . 

7C XL ( IA ) = P ITCH 

CO 80 JA = MXBI.MX6U 

NU(IA) = INT ( {Xl!( IA )+DTLR)/HB) 

IF! XU I IA > .GT .-DTLR ) NU ( I A ) =NL < I A) + 1 
NL!IA) = 1NTIIXLI IA)-OILR)/HB) 

80 1F!XL< IA) .LT.DTLR) NLCIA) = NLUAI-I 


CO 9C IA = MXBOP 1 *MX 
NU( IA ) = NBUO ' 

NL 1 I A ) = NBBO 
XU! IA ) = STGk 


9C 

300 

998 
1000 
ICO 5 
1010 
1020 
HOC 

1110 

1120 


XL! IA ) = 
RETURN 
CALL CEUEKR 
RETURN 
FORMA 1 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 


P ITCH + STGR 


146H MX,N8BI , NU S P , N L SP , NR S P , CR N 1 NT IS TOC LARGE*) 


< IF1) 

! 1615 ) 

< IX, 1617 ) 

( 8F10-5 ) 

( IX, 8016. 7) 

( 7X, 3FGAM, 14X, 2HAR ,13X,3HTI P ,12 X ,5 H RHO I P , 12X ,4HWTFL,11X, 


1130 
1140 
1 15C 
1 160 
1 1 7C 
1180 
1 190 
120C 
121C 


1 5 HO MEGA, 13X, 1HW) 

FuRMAT ! 6 X , 5 HC HQRO ,12X,4HSTGR,11X,5HBETAI ,1 1 X ,5 HeET AO ) 

FORMAT (8X,2hRI, 13X.4HALUI ,12X,4HALLI ,13 X ,2HR0 , 13X , 4HALU0 , 12X , 
1 4h ALLO ) 

N8BI NUSF NLSP NRSP NBL NINT) 


FO RM AT 
FOkM A T 
FORMAT 
FORMAT 
FORM A T 
FORMAT 
FORMAT 
FORMAT 
FuRMAT 
1EL ) 

1220 FORMAT 
1 BETA 
20 *W 


1 4 £H MXB I M XB 0 MX 
< 1CX, 6HML ARRAY) 

( 1CX, 1CHXSPU ARRAY) 

( 1CX, EHML ARRA Y) 

( 1CX, 1CHXSPL ARRAY) 
( 1CX, 8HMR ARRA Y ) 

( 1CX, 1CHKMSP ARRAY) 
l 1CX, 1CFCESP ARRA Y) 
(57H BLDATA NULAKI 


ERPRT STRFN SLCRD ARPRT INTV EL SURV 


! 5E1LAMBDA = , G 14 . 6 ,/ 12 X , 6 8HFRE E STRE AM 
CORRECTED TO BLADE 
BLADE LE OR TE 


MAXIMUM 


VALUE 

CRI T1 CA L/l 3X,63HVELCCIT Y FOR RH 

VELCCITY/8H INLET .4G18.7/8H OUT 


3LET , AG18.7) 

1230 FORMAT ( 1 FL , I0X , 2 8HCA LC OLA TEU PRUGRAM CONSTANTS/ 7X ,2HHA , 14X , 2HHB, 

1 14X , 2H FU , 14X , 2HhL , 13X , 5HP1 TCH ,12 X ,4HNBBt ,6 X,4HNBU0/ 1X,5G16.7, 

2 2110 

1240 FORMAT (78hL CALTICi', - H8*RM ( MXBC) LESS THAN RU LESS THAN HA MAY N 
10T GIVE CORRECT RESULTS) 

1250 FORMAT ( 1 FI , 13X , 48HS TrE AM SHEET COORDINATES AND THICKNESS TABLE / 

1 7X, 1FM, 14X,1HR,13X,3HSAL,13X,1HB,12X,5HD8/0M / (5G15.5)) 

1260 FORMAT { IF 1 » 13X, 27HBLAUE DATA AT SPLINE POINTS /18 X ,16 HUPPER SU 
1RFACE ) 

1265 FORMAT { 18X, 16HL0WER SURFACE) 
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1270 FORMAT ( 7X , IBM , 10X , 5HTHE TA , 1 CX .10HDER I VATI VE ,5 X ,1 0H2N0 OERIV. / 
1 (4G15.5) ) 

128C FORMAT ( 1 F 1 , 13X , 22HBLADE COORDINATE TABLE/ 7 X ,1 HM , 14X ,2HXU , 11X , 

1 5HUXC2U, 12X, 2HXL , 11 X, 5HD XD2L / (5G15.5) ) 

140C FORMAT (22FL WEIGHT FLOW REDUCED T0.G14.6.7H KG/SEC) 

END 


SUBROUTINE COEP(NXNI) 

REAL K »K1»K2»K3,K4,K5, LMA X , LPI N , LAMBDA .MU.PL.MR.PXU ,MXL,MPL 
INTEGER SRW, FIR ST, CASE ,BLUA TA ,ERPR I , S T RF N ,S LCRD ,SLPLT , ARPRT ,SURV EL 
COMMON SRW, ITER ,GAM,AR ,T I P.RHOI P , WlTF L , GREG A , L AMBDA , CP , EX PON, P IT CH, 

1 CHORC, STGR»8ETAI,BETA0,UTLR,RI , A LUI »ALLi , RC , A LUC , ALLC , 

2 MXB 1,MXB0 , M X , N8B I , NUSP ,NLSP ,NRSP » M NT , VTCL , 

3 BL DATA, NO LA KI,ERPRT,STRFN,S LCRD, ARPRT ,1 MVEL .SURVEL, 

4 MO ( 20) , XSP L ( 50 ) , ML (50, XSPLI50) , MR ( SO ) , RMS P 150 ) , 8ES P 1 50 ) , 

5 w, mK, TOLER, BOA ,BD0,U(250G) ,A(2500,4) ,K(25GG) , FH0(250G) , 

6 CXO/UI ICG I.DXDZLIICC) .SLOPE (50) ,EMC (50) ,SLLPE (50 ) ,EML(50), 

7 RM( ICO ),BE( ICC) , SAL ( IOC) ,Xt ( 100) , X L (100 ) , RMU (100 ) , RML ( 100 ) , 

8 NU( ICO), NL( ICC), Cl NT (11) , XI NT (11) .PPL (100) , 

9 hA, FB.NXN , MXB IM1.MXB0P 1 , JU.JL.HU, HL.NBBC.NBUC ,NCH, 

1 I8TE1, IBTE2, ITE2.HA1TE.HAME ,CNI ,LNC,TWw»lTERA, 

2 RHO U( ICC) ,RHUL< ICC ) .RHOCT (ICC) , RHC L T ( 1 00 ) , BEU ( 100 ) ,BEL( 100 ), 

3 AAA (IOC) 

C THE FOLLOWING VARIABLES ARE ALL IN CCPMCN BY EQUIVALENCE 

DIMENSION WM(25CC),WX(250C) ,V(2500) , BE TA (2500 ) ,S L ( 1 100 ) 
EQUIVALENCE ( A ( 1 , 1 ) , WM { 1 ) ) , ( A ( 1 ,2 ) , WX ( 1 ) ) , ( A ( 1 ,3 ) , V ( 1 ) ) , 

1 ( A ( 1.4). BETA( 1)),(K(14G1) , S L ( 1 ) ) 

ITER = ITER + 1 

IF( GAM .EQ . 1 . 5. AND .AR.EQ. 100C..AND. TI P. EC.l . E6 ) GC TC 94 
IF( ITEK.N E. 1 .AND. I TER .NE .2) GO TC 92 
NULAK I = NULAK I- 1 
ERPRT = ERPRT-1 
STRFN = STRFN-1 
SLCRC = SLCRO-1 
ARPRT = ARPR T- 1 
INTVEL = INTVEL-1 
SURVEL = SURVEL-1 
92 IF( ITER.NE.0) GO TO 94 
NULAK I = NULAK 1 + 2 
ERPRT = ERPRT+2 
STRFN = STRFN+2 
SLCRG = SLCRU + 2 
ARPRT = ARPRT+2 
INTVEL = INTVEL+2 
SURVEL = SURVEL + 2 
94 I = C 
INL = C 
INU = C 

IF(BLCATA.GT.C) WRI TE (6, 144G) 
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c 

C CALCULATE COEFFICIENTS 
C 

CO 20 C IA = 1 » M X 
XI = FLOATINLI I A ) )*HB 
I F( 1A.C-T.MXBO) XI = STGR 
NBC = NLC 1A >-NL( IA )+l 
K 5 = C . 

IBUU = N U ( IA— 1 I— NO ( I A ) 

I60L = NL ( IA )— NL ( IA + 1) 

CO 20 C 18=1, NBB 
I = I + I 

1 Ft I.C-T .2500 ) GO TO 200 
1F( IA .NE. 1) GO TO ICC 
A( I, 1 ) = C . 

At 1,21 = C. 

A( I, 3 ) = C . 

A( i, 4 ) = 1 . 

K( I ) = FA*UN I 
GU TO ISC 

IOC IF(IA.NE.MX) GO TO 11G 
At I, 1) = C. 

A < 1,2) = C . 

A( 1, 4 ) = C . 

A < 1,3) = 1 . 

KID = FA * UNO 
GO TO ISO 
11C K 1 = 0 . 

K 2 = C • 

K3 - C . 

K 4 = C « 

H 3 = FA 
Fi4 = FA 

I Ft (NU( IA )+IB ) .LE.NUt IA-III CALL 1 N TP Lit XU , I A , XI ,H3 ,0 , K3 , MU , XS PU , 

1 ALU I, ALOO ,NCSP , SLUPE , EMU ,RI ,RO, CHORD , STGR, PITCH , l.,HA,hB, 

2 MXB I,M XBG »RM l MXB I ) ,RM(MXBC) ) 

IFt tNUUA ) + IB ) .LE .NUl I A+ 1)1 CALL I NT P L I ( XU , I A , X I ,H4 , 1 , K4 , MU , XS PU , 

1 ALU 1, ALUO ,NL SP , SLOPE ,EMU ,RI ,RG , CHORD, STGR, PITCH , l.,HA,hB» 

2 MXB l,MXBO,RMIMXttI ) ,RM(MXB0) ) 

IFttNUtIA ) + IB) .GT.tNL(lA-l)Fl ) ) CALL INTPLl(XL,IA,Xi,H3,0,K3,KL, 

1 XSPL, ALL I,ALLO*NLSP»SLLPE ,E ML , R I , RC, CHORD, STGR, PITCH, -1., HA, HB, 

2 MXB I, MX BO ,RM(MXBI ) ,RMIMXBC) ) 

IFt ( N U ( 1 A ) + 16 1 . GT . ( NL ( IA + 1 ) + 1 ) ) CALL I NTP LI (XL, I A, XI ,H4, 1 , K4, ML, 

1 XSPL, ALL I ,ALLO, NLSP , SLLPE ,E ML »RI , RC, CHORD, STGR, PITCH, -1., HA, HB, 

2 MXB 1, MX BO, RMtMXB I 1 , RM ( M XB C 1 ) 

IFt IA-MXBO+1 » 140 , 12C , 13C 

C SPECIAL CALCULATION OF COEFFICIENTS AT TRAILING EDGE LF LOWER SURFACE OF 
C UPPER BLADE 

120 IFtXl .LT .P ITCH+STGR— RO/RMtMXBO) ) GC TC 140 
IFtXL.CT.P 1TCH + STGR) GO TO 140 
IFt IB+NUt IA )— l.GT.NBBO) GO TC 14C 
HAMRO = HA-RO 

XL ( MX BO I = P ITCH + STGR-RO/RM tMXBG) 

CALL I NTP Lit XL,IA,X1,H4,1,K4,ML,XSPL,ALLI , ALLO.NLSP, SLLPE, EML, 

1 R I, RO, CHORD, STGR ,PI TCH ,-l . ,H A MRC , H B , MXBI , MXBO , RM t MX B I ) , RMt MX BC ) ) 
HANTE = H4 
IBTE2 = N L ( IA J + IB-l 
XL t MX BO ) = P ITCh+STGR 
IBLL = IBOL+1 
GO TO 140 
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13C I F< IA .C-T.MXBO) £0 TO 140 

IFT IBTE2.EU. 10CC ) GO TO 14C 

1FTX1.LT.PITCH+STGR-R0/RMTMXB0) ) GO TC 140 
B = XLTMXBO-J) 

XL ( MX BO— 1 I = P ITCH+STGR-RO/RMTMXBC) 

CALL INTPL1T XL , I A , X1,H3,0,K3 ,ML ,XSP L.ALLI , ALLC , NLS P ,S LLPE, EML, 

1 R I » RO » CHORD » STGRiPI TCH »— 1. tRC,HB,RXBI , MXBO , RR T MXBI ) .RMTRXBO) ) 
HA1TE = h 3 
IBTE1 * NL( I A ) + IB- 1 
ITE2 = I 
XL T MX8G- 1 ) = B 
140 IFT IB.EQ. I > K 1 = 1. 

I Ft IB.EQ .NBB I K2 = 1. 

Hi - HE*RM( IA ) 

h 2 = HI 

II = 1-1 

12 = HI 

13 = I-NL 1 IA- 1 ) +N L T I A )- 1 

14 = HNLT IA 1-NUIA + ll + l 
IF! IA .GE.MXB I) GO TO 142 
IF( IB. tU. 1) 11 = 1 1+NBB 

I F C IB . EQ .N8B ) 12 = I2-NBB 

GO TU 148 

142 I Ft IA .CT .MXBO ) GU TO 144 

IF! IB.EQ. 1) HI = (Xl-XUTIA) )*RM(IA) 

I FI IB . Ey .NBB ) H2 = ( XL ( I A ) - Xl ) *RM ( I A ) 

GU TO 148 

144 IF( IB .NE. 1 ) GO TO 146 

II = Il+NBB 

hi = HL*RM T I A ) 

H2 = hL*RMT IA ) 

GU TU 148 

146 IFT IB.EQ. 2) HI = HU*RM(IA) 

I FT I B .RE .NBB 1 GO TO 148 
12 = I2-NeB 
H2 = hL*RM( IA ) 

148 RHQ 1 = RHO ( I 1) 

Rhu 2 = RHO ( 12 1 
RHU 3 = RhO( 13 ) 

RHU 4 = R hO ( 14) 

BE3 = BET IA- 1 ) 

BE4 = BET I A ♦ 1 ) 

IFT T IA .GE .MXB I .AND. IA .LE .MXBO) . AND. I B.EC. 1 ) RHG1 = RHOUTIA) 

IFT ( I A .GE .MXB I .AND .IA .LE.MXBG) . AND. I B. EG. NBB) RHC2 = RHCLTIA) 

IFTKS.LT. .5) GO TO 16C 

IFTK3.LT . .5) GO TO 15C 

INL = 1NL+1 

BE 3 = BEL ( INL ) 

RHO 3 = RHOLTT INL ) 

I FT B E 3 .NE .C. ) GO TU ISC 
B = MPLT IA )-F2 

CALL SPLINT (MR ,BESP ,NRSP ,B ,1 ,BE 1(1 NL) ,AAA) 

B E 3 = BEL ( INL ) 

IFT BLCATA.GT.C ) HR I TE ( 6 , 1 46 C ) B » I NL 
15C IFTK4.LT. .£) GU TO 18C 
INL = INL +1 

I TP = INL — IBDL + 1 + 2* I NB B- IB) 

BE4 = eEL ( I TP ) 

RHO 4 = RHULTT IIP ) 

IFTBE4.NE.G.) GU TO 16C 
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b = M P L ( IA } + F4 

CALL SPLINT ( MR , BE SP , NR SP , B , 1 ,BE L ( i TP) , A A A) 

BE4 = e£L { I TP ) 

IF( BLCATA .GT .0 ) WR I TE (6*1460 B.INL.ITP 
GU TU 180 
16C CUNTINLE 

IF( K 3 • LT . .5) GU TU 17C 
1NU = 1NU + 1 
ITP = 1NU+ IBDL + 1-2* 18 
6E3 = EEC ( ITP ) 

RHU3 = RHJLT ( ITP ) 

IF( BE2 .NE.C. > GO TO 17G 
B = MPL( I A )— H 3 

CALL SPLINT ( MR , BESP ,NR SP , B , 1 ,BEL( I TP) ,AAAI 
BE3 = EEUI ITP J 

I Ft BLCATA .GT .C ] wR I TE ( 6, 1450) B,INU,ITP 
17C IFIK4.LT. .5) GU TU I8G 
INU = 1NU+1 

BE4 = EEUI INC) 

RhU 4 = RHUUTI INC ) 

I F ( BE4 .NE .C. ) GU TU 18C 
B = MPL( I A ) + H4 

CALL SPLINT (MR, BESP, NRSP,B,1, BEHIND , AAA) 

BE4 = BEU ( INU ) 

I Ft BLCATA .GT .C ) WK l TE (6,1450 B,INU 
18C CUNTINLE 

A 12 = 2./F1/H2 
A34 = 2./F3/H4 
AH = A 12 -» A 3 4 

B 12 = (RHQ 2— KHU 1 ) /RHO (I) /(HI +H 21 

634 = (BE4*RhU4— BE 3*RHU3 ) /BE 1 1 A ) /RFIO (I ) / (H3+H4) - SAL(IA)/RP( IA) 
At 1,11 = ( 2./H1+B 12)/IH1+H2) /At 1 
A( 1,2) = A12/A1I— A(I ,1) 

At 1,3) = ( 2./H 2+B34)/IH3+H4) /At I 
A( 1,4) = A34/A I I-At I ,3) 

K ( I ) = — T W W*BE ( IA )*RHU ( I ) * SAL ( I A) /All 
IFIK3.LT. .5. AN0.K4.LT. .5) K5 = 1. 

I Ft IA.LT.MXB l.UR.IA.GT.MXBU) GC TU 185 

KID = K( I ) + K5*(K2*A ( I ,2)*K3*A (I ,3) + K4*A (1 ,4) ) 

IFIKi .GT. .5) A( I ,1) = C. 

IFIK2 .C-T . .5 ) A( I ,2 ) = C. 

18 £ IFIK3.CT..5) A<I, 3) = C. 

IFIK4.CT. .5) A(I ,4) = C. 

19 C XI = FLOATING! IAJ+IB)*HB 

IF (1A.GE.MXBI .AND .IA .LE.MXBC) GO 1C 200 
IF(IB.EU.l) K(I) =-A(I,l) + MI) 

IF( IB .EQ.NBB ) K(I) = A(I,2)+K(I) 

20 C CUNTINLE 

I F ( IT ER.EQ. 1 )RR I TE 16,1470 INU,INL 
IF( INU .GT . 100. OK. INL.GT. ICO) GG TO 310 

1FUITER .GT. D.AND.dNU.NE. JL.CR.1M.NE. JL> ) GC TC 320 
NXN = 1 

WRITE (6, 1420 NXN 

IFINXN .GT.25CC) WR I TE ( 6 , 1430 ) 

NXN 1 = NXN 

1F( BLCATA .LE .C ) GO TU 21G 

WRITE (6,1410 (NU(lA),NL(lA) ,IA = 1,FX> 

21G BLCATA = C 

RETURN 


non 


310 CALL CEOERR OSH INU AND INL MOST BE LESS THAN 100$) 

220 CALL 0E0ERR ( 44H INC NOT EQUAL TO JO CR INL NOT EQUAL TC JL$) 
RETURN 

1410 FORMAT ( 16 HL LIST OF NU AND NL / (2110) 

1420 FORMAT 1 lhL, 5X , 32HNUMBER OF INTERIOR MESH POINTS = ,I5) 

1430 FORMAT 1 76HKT HE NUMBER OF UNKNOMN MESH PCI NTS EXCEEDS 2500, A COARS 
1ER GR IC MOST BE USED) 

1440 FORMAT OlhL M INU INL I TP) 

1450 FORMAT ( IX , G13 .4 , 1 5 , 1 10 ) 

146C FORMAT ( IX , G 13 .4 , 1 1C , I 5 ) 

1470 FORMAT ( 6HL IN U =, 13, 5X.5HINL =,I3) 

END 


SUBROUTINE SOR 

REAL K,K1,K2,K3»K4»K5» LMAX , LMI N , LAMBDA , HU , ML , MR , MXU ,MXL,MPL 
INTEGER SRW, FIRST, CASE ,BLO A TA ,E RPRT , ST RFN ,SLC RD ,SLPLT , ARPRT .SURVEL 
COMMON SRH, ITER , GAM, AR ,TIP,RHOI P,WTFL , CMEG A , L AM eDA , CP , EXPON, PIT CH, 

1 CHORD, STGR,BETAl,BETAO,OTLR,RI ,ALUI , ALLI ,RC , ALUC , ALLC , 

2 MX B I,M XBO , MX , NBB I , NU SP , NL SP , NR SP , NI NT , VT UL , 

3 BL DATA , NULAKI , ERPR T , STRF N , SLCRD , AR PRT , I NTVEL ,SURVEL, 

4 MU( SO) ,XSPU( 50) ,ML I 50) ,XSPL(50) ,MR(50) ,RMSP(50) ,BESPI50), 

£ H, MR, TOLER ,BDA ,BDO ,0(2500) »A (2500 ,4) ,K(2500 ) ,RH0(2500 ) , 

6 OXOZUI 100 ),DXDZL ( 100) ,SLUPE (50) ,EMU(50) ,SLLPE(50) ,EML(50), 

7 RM( 10C),BE( 100), SAL ( 100 ) ,X0 ( 100) ,XL(100) , RMU (100 ) , RML ( 100 ) , 

8 NU( ICO ) ,NL ( 1C0),UINT(11) ,XINT(11) ,MPL(100) , 

g HA, EB, NXN ,MXBIM1,MXB0P1 , JU, JL ,HU ,HL , NBBC ,NBUC ,NCH , 

1 IBTE1, IBTE2, ITE2,HA1TE ,HANTE ,UNI ,LN0 ,TWM,ITEFA, 

2 RHOUI ICC) ,R HOL ( 100) ,RHOUT(100) , RHCLT (100) , BEU (100) ,BEL( 100 ), 

2 AAA (IOC) 

C THE FOLLOWING VARIABLES ARE ALL IN CCMMCN BY EQUIVALENCE 

DIMENSION MM (2500), WX( 25 00) , V ( 2 500) , BETA (2500 ) ,SL ( 1100 ) 
EQUIVALENCE (A(1,1),MM(1)),(A(1,2),MX(1)),(A(1,3),V(1)), 

1 ( A ( 1,4), 8ETA( 1) ) , (M14C1) ,SL(1) ) 

DIMENSION IU( ICC ),IL( ICO) 

EQUIVALENCE (K ( 1 10 1 ) , I U( 1 ) ) , ( K (1 201 ) , I L (1 ) ) 

ESTIMATE MB E ST 

NUTEMP = NULAKI 
IF(M.GE.l.) GO TO 225 
19C DO 20C 1=1, NXN 
20 C U( I ) = 1. 

MM AX = 2. 

21C I = C 

MM AX 1 = WMAX 
LMAX = 0. 

LMIN = 1. 

00 22C IA = 1,MX 
NBB = NL( IA )-NU( IA ) + l 
CO 22C IB =1 ,NBB 

1 = i*: 
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n n n non 


11 = 1-1 

12 = m 

IF( ( ( 1A.LT.MXB I ).0K. ( IA.GT.MXBL) ) .AND. < I B. EQ. 1) 1 11= Il+NBB 
I Ft ( ( IA.LT .MXB I ).UR.(IA.GT.MXBO) ). AND. (IB. Eg. NBB) ) 12= I2-NBB 

13 = I-NL( 1A-1)+NLUA)-1 

14 = 1 +Nl ( IA )-NU( IA + 1 ) + l 

UNfcW = A( I, i)*U!(Il)+A(i ,2) *U(I 2)+A (I ,3 > *U ( 1 3 ) + A ( I ,4 ) *U ( 14 ) 

HAT IU = UNEW/LId) 

LM AX = AMAX1(RAT10,LMAX) 

LMIN = AM INl(RATiO.LMIN) 

215 U< I) = UN EH 
220 CONTINUE 

WM AX = 2./( 1.+SQR1 (ABSU.-LMAX) )l 
WM IN = 2./( l. + SQRTd.-LMIN) ) 

WRITE (6,150 0 WMAX,WMIN,LMAX,LMIN 

I F( ( ( WMAX 1-WMAX ) .GT. WR ) .OR. ( WMAX.G I. (2.-100.*WR) ) ) GO TC 210 
W = WM AX 

CALCULATE INITIAL SOLUTION ESTIMATE 
225 I = C 

IF( ITER.NE.l) GO TO 26C 
DU 23 C IA = 1 , M XB IM 1 
NBB = NL( IA )-N U ( I A ) +1 
CO 23 C IB = 1 » NBB 

1 = 1+] 

U( 11 = FLOA T ( IB— 1 ) /F LOA T ( NBB ) 

23C CONTINUE 

00 24 C l A =MXB I ,MXBO 
NBB = NL ( IA J-NUl IA J+i 
CO 24 C IB = 1 , NBB 

1 = 1 +] 

J = N U ( IA J * I B— 1 

UII) = ( HB*FLOAT( J I-XU(IA) ) / ( XL ( I A ) -XU ( I A I) 

240 CONTINUE 

CO 250 IA =MXBUP 1 , M X 
NBB = NL ( 1 A ) — N U ( IA ) + l 
CU 2 5 C 18 = 1, NBB 
1 = 1+1 

Util = FLOA T ( IB- 1 ) /FLOA T (NBB ) 

25C CONTINUE 

SCUVE MATRIX EQUATION BY SUR 
260 I = C 

IFtNUTEMR .GT.C) WRITE (6, 1450) 

ERROR = 0 . 

CO 27 C IA=1,MX 
NBB = NL( IA )-NU( IA ) + l 
CO 2 7 C IB = 1 , NBB 
1 = 1 +] 

II = 1-1 
12 = 1+1 

1F( ( ( IA.L T .MXB I ) .OR. ( IA.GT.MXBC) ) .AND. (IB. EQ. 1) ) 11 = ii + NBB 

IFI (( IA.LT .MXBI ) .OR .( IA.GT.MXBO) ) .AND. (IB. EC. NBB) ) 12= 12-NBB 

13 = 1-NL ( 1A-1 ) +NU( IA )— 1 

14 = I +NL ( IA )— NU( IA + ll + l 

CHANGE = W*(M I )— U( I ) +A ( I , 1)*U(I1)+-A (1 ,2) *0(12) +A(1 ,3) *U ( 13) + 
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1 A( 1,4 )*U( 14)) 

ERROR = AM AX 1 ( ERROR ,AB S (C HANbE ) ) 

U( 1 ) = U( I )+C RANGE 
IF(NUT EMP .LE .G ) GO TO 27C 
IF(lA.EU.l) 1 3 = C 
I F( IA.EQ.MX) I 4 = C 

WRITE (6, 1460) IA , I , (A ( I , J ) , J=1 , 4 ) ,1 1 , 12 , 1 3 , 1 4 , K ( I ) 
27C CONTINUE 
NUTEMP = C 

IF( ERPRT.CT.O) WRITE (6,1510) ERROR 
IF( ERROR. GT. TOLER ) GO TO 260 
IF( STKFN.GT.C ) NR I TE ( 6 , 152 C ) 

LAST = C 
DO 26C IA=i,MX 


IF( STRFN.GT.C ) NR I TE (6,1525) IA 

FIRST = LAST ♦ 1 

LAST = F1RST+NL(IA)— NU(IA) 

IU( IA )=FIRST 


A( 1,4) 


1L( 1 A ) =LA S T 

26 0 iFl STRfN.GT.C ) WRITE (6,1530 ( L (I ) , I =F I PST , LAST ) 

RETURN 

1450 FORMAT (8EF1 IA I A(I,1) A(I,2) A(1, 3) 

111 12 13 14 K ( I > ) 

1 46 C FURMAT ( 2 X , I 3 , 1 6 , 4F 10 . 5 ,4 1 7 ,F 1 C . 5) 

1500 FORMAT ( 7h WMAX = ,F 9 . 6 , 5X , 6H WMI N = ,F 9 . 6 ,5 X ,6H LMAX = ,F9.6,5X, 
1 6HLMIN = , F 9 . 6 ) 

1 5 1C FORMAT ( EE ERROR =,F11.6) 

1520 FURMAT ( IF 1, 10X , 22HSTREAM FUNCTION VALUES) 

1525 FORMAT (5F IA = , 1 3 ) 

1530 FORMAT (2X,1CF13.£) 

ENC 


SUBROUTINE SLAXVL 

REAL K ,K1,K2»K3,K4,K5,LMAX,LRIN, LAMBDA ,MU , ML , MR , MXU ,MXL,MPL 
IN TtGER SR W, FIR ST, CASE ,B LU A TA , E RPR T , S T RF N , SLCRO , SLPLT , ARPRT ,SU RV EL 
COMMON SR W, ITER ,GAM , AR , T 1 P ,RHOl P ,W IF L , C MEG A , L AM EC A ,CP , EXPCN, P IT CH, 

1 CHORD, S TOR, BETA I ,BE TAO ,0 T LR ,RI , A LLI , AL LI , RC , A LUO , ALLC, 

2 MXB I, MX BO ,MX,NBBI , NUSP , NL SP , NR SP , M NT , VTCL , 

3 BLCATA , NULAKI ,ERPRT , STRF N ,SLC RD , ARPRT , I NTVEL ,SURVEL, 

4 MU( 50) ,XSPL( 50) ,ML ( 50) ,XSPL(50) , M R (50) , RMS P ( 50 ) , 8 ES P (50 ) , 

5 W, MR, TOLER, BOA ,BDD, 0(2500 ,A(2500,4) ,*(2500) ,RHC(2500), 

6 OX 02U ( ICO ),DX02L ( 1UG) ,SLUPE (5 0) »E ML (50) ,SLLFE (50 ) , EMU 50), 

7 RM( ICO ),BE( ICC) , SA L ( IOC) , XL ( 100) , X L l 1 00 ) , RMU (100 ) , RML ( 100 ) , 

8 NU( 1C0),NL( 1CC),UINT(U) , XI NT (11) ,MPL(100) , 

9 HA, FB,NXN,MXBIM1,MXB0P1 , JU , JL ,H U ,H L , N8 BC , NBUC ,NCH, 

1 IoT El, IBTE2, I TE2,HA1TE .HANTE ,LNI , L NC , T n W , I T E R A , 

2 KHO U( ICC) ,R HUL ( 100) ,RH00T(1CG) ,RHCLT(100) ,BEU(100) , BEL (100), 

3 AAA (IOC) 

C TEE FULL CW ING VARIABLES ARE ALL IN CCMMLN BV EQUIVALENCE 

DIMENSION WM(25GU),WX(25CC)»V(2500) , BE TA (2500 ) ,SL ( i 100 ) 

OIMENS IUN WMU( ICC) , WML( ICO) ,WXL (100) ,WXL( 100) , MXU (100 ) , MX L ( 100 ), 
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o o o r> n n 


1 BETAU11CC), DETALt ICC) ,WU(ICC) .hLtlGO) , 

2 USP(ICC), ICt ICOfILtICG) 

EQUIVALENCE I A I 1 , 1 ) , WM 1 1 ) ) , I A 1 1 ,2 ) » in X ( 1 ) I , I A 1 1 »3 ) , V 1 1 ) ) , 

1 IA1],4), BETA) 1)),IKU*#C1) tSL(l) ) 

EQUIVALENCE (K(l) ,WMUI 1) ) , I K (1 Cl ) , WML 1 1 ) ) , ( K( 20 1 ) ,WXUl 1 ) ) , 

1 (K130I ),WXL( 1)),<K14C1) ,MXU (II) , 1X1501) ,MXLI1)) , 

2 ( X ( 6C1 ) » BETAUI i) ) ,(M7C1) , BETAL II )) , I K (80 1 ) , WU 1 1 ) ) , 

3 (K< SCI ),WL( i ) ) , (M 1CC1) ,LSP(1) ) , IK (1101 ) ,IU(1) I , 

4 < K f 1201 )> IL(l) ) 

DIMENS ION XBBI 52) ,WMSPI52) ,KKK<24) ,P(I1) 

DATA IKKM J),J=4,24,2)/11*1H*/ 

CALCULATE STREAMLINE LUCAT10N — UPSTREAM 


11=1 
I = 0 

0 EL IN T = 1 ./FLOAT (NINT) 

NBB = NL ( 1 )+l 

XBBI 1 ) = C. 

IF( SLCRD.GT.C) WR I TE ( 6 • 1 55C ) 

CU 320 IB = 1 »NBB 
32C XBBI IB +1 ) = XbB I I B ) +HB 
CU 35C IA = l.MXBIMl 

U INTI 1 ) = AINT10II+1)/UELINT)*DEL1NT 
IF! U< I-»l) .GT.C.) U I NT 1 1 ) = LINT (1 ) + DE LI NT 
CU 33 C JB = 2,N INT 

33C UINT1JE) = UINTI JB-ll+DELINT 
CU 34 C IB = 1 , NB B 

1 = 1 +] 

34C USPI IB ) = 01 1 ) 

NSP = NBB+1 
USPINSP) = USPI 11 + 1. 

IF! SLCRC.C-T.C) CALL SPLINT I USP ,XBB , N S P , 01 NT , NI NT , X I NT , AAA ) 

CALL SPLINE l X6B , USP ,N SP , WM 1 1 1 ) » AAA) 

I 1= U+NL1 IA ) +1 

IF! SLCRD.GT.C) WK I TE I 6 « 1 f>4C ) MPLII A ) , I 01 NT I J) , XI NT I J ) , J= 1 , N INT ) 
CU 345 JB=1,NINT 
J = MX + IJB-1 ) + I A 
SLI J ) = X INT I JB ) 

245 CUNT IN LE 

J 1 = MX*N INT + I A 
SHJ1) = SLIJ) 

25C CONTINUE 

CALCULATE SIREAMLINE LUCATICN — BLADE 

J U =MX B 1- 1 
NINT = NINT+1 
UINTI 1 ) = C. 

CO 36 C JB = 2»NINT 
36C UINTIJE) = 0 INTI JB-U+DELI NT 
USPI 1 ) = C . 

DU 38 C I A =MXB I *M XBQ 
XBBI 1 ) = XUI IA ) 

NBB = NL t IA ) — N L I IA) + 1 
XBBI2) = FLOATINUI IA) )*HB 
CU 37 C IB = 1 > NBB 
1 = 1+1 

USPI IB+1) = Oil) 
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37C XBBIIB+l) = FLOAT ( IB- 1 )*HB + XBB (2) 

NSP = NBB+2 
USP(NSP) = 1. 

XBBlNSP) = XL ( IA ) 

I FC SLCRD.GT.C) CALL SPLINT { LlSP ,XBB ,NSP , LI NT . N1 NT ,X INT , AAA ) 

CALL SPL IN£{ XBB ,USP ,NSP ,WMSP ,AAA) 

DO 375 IB = It NBB 
WM ( II )=WMSP( IB ♦ 1 ) 

275 11=11+1 
JU=JU + ] 

WMUIJU > = WMSPI 1) 

MMLC J U ) = WMSPINSP) 

IF! SLCRD.GT.C) WR I T£ 1 6 , 1 54C ) MPL (I A! , (LI NT( J) »XINT I J1 , J=1,NINT ) 
DO 38C JB = It NIN T 
J = MX*(JB-1) + IA 
SLIJ ) = X INT ( JB ) 

28C CONT IN LE 
C 

C CALCULATE STREAMLINE LOCATION -- DOWNSTREAM 
C 

NINT = NINT-1 

NBB = NLIMXBUP 1 )-NU( M XBGP1 ) + 1 
NSP = NBB +1 
X BB{ 1 ) = STGR 
XBBI2) = STGR+HO 
CO 39C 18=3, NBB 
290 XBB(IB) = XBBIIB-1I+HB 
XB6( NSP I = STGR+PITCH 
DO 42C IA =MXBOP 1 ,MX 

UINT(l) = AINTIUII+l) /DELI NTI *DELI NT 
IF( U( I +1 ) .GT .0 . ) UINT(l) = UINTCD+OELINT 
DO 40 C JB=2»NINT 

40C UINTIJE) = UINTI JB-U+DELINT 
DO 41C IB = 1, NBB 
I = I+I 

410 USP< IB ) = U( I) 

USPINSP) = USPID + 1. 

I FC SLCRD.GT.O) CALL SPLINT ( USP ,XBB , N SP , LI NT , NI NT , X I NT , AAA I 
CALL SPLINE! XBB «USP,NSP,WM (ID , AAA ) 

1 1= I1+NL( IA )— NL( IA 1 + 1 

IF! SLCRD.GT.C) WR I TE I 6 » 1540 MPL (I A) , (III NT (J) , X I NT ( J ) , J= 1 , N I NT ) 
DO 415 JB = I* N IN T 
J = MX*(JB-1)+IA 
SLIJ) = X INT ( JB ) 

415 CONTINUE 

J1 = MX*N INT + IA 
SL(J1) = SLIJ) 

42C CONTINUE 
C 

C PLOT STREAMLINES 
C 

IF! SLCRC.LE.C) GO TO 460 
C CALCULATE PLOTTING PARAMETERS 
2MIN = MPL I 1 ) 

XMAX = XL! 1) 

XM IN = XU I 1) 

DO 430 IA = 2, MX 

XMAX = AM AX II XMAX, XL I IA ) ) 
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430 XMIN = AM INK XM1N,XU(IA) ) 

OX = XMAX-XMIN 
X FACT = 2 . 

44C IF{ OX .CT. 10. ) GO TO 45C 
OX s OX* 1C. 

X FACT = X FAC T + 1 . 

GO TO 440 

45C IF( CX.LE . ICO. > GO TO 460 
OX = OX/1C. 

XFACT = XFACT-1. 

GO TO 450 

46 C OX = A 1NT ( DX + 1 . ) 

DZ = A 1NT< 5.*DX/3.*RM(MXBI 1 ) 

Z FACT = XFACT 

ZM IN = AINT( ZMIN*10.**ZFACT) 

XM IN = AINT< XMIN*1G.**XFACT) 

KKK ( 1 ) = 45 
KKM 2) = NINT+1 
P( 1) = 1. 

P( 5 1 = C. 

P( 6) = 6.-ZFACT 
P< 7 1 = ZM IN 
P( 8 ) = OZ 
P( 9) = 6. -XFACT 
P< 10) = XM IN 

P< 11) = CX 
KKM 3 ) = MX 
WRITEt 6, 1530 ) 

CALL PLGTMY ( MPL , SL ,KKK , P ) 

WR ITE (6, 1 56 C 1 
CG 47C IA = 1 » M X 

47 C MPL ( I A 1 = FLGAT(IA-MXBI)*HA 
48C LAST = C 

CO 482 IA = 1 » M X 
FIRST = LAST + 1 
LAST = FIRST +NL(IA)—NU(IA) 

CQ 482 I = FIRST,LAST 
482 WM(I) = WM ( 1 )*WTFL/BE(1A >/RM(IA) 

OU 484 l A =MXB I ,M XBU 

wMUtiAl = WML( IA 1* WTFL/BE (I A) /RM(I A ) 

484 WML ( 1 A ) = HML< I A )* WTF L /BE ( I A ) /R M < I A > 

I Ft ARPPT .LE .C ) KETLRN 
WRITE (6, 1600 
L AST = C 

CG 49 C IA = 1, MX 
FIRST =LAST + 1 

LAST = FIRST+NL( IAJ-NU(IA) 

49C WRITE (6,1020 ( WM ill , I =F I R S T , L A ST) 

WRITE (6, 1 6 1 C 1 

WRITE (6, 1C2C) IWMUIA), IA=MXBI ,MXBC> 

WRITE (6, 1620 

WRITE (6, 1C2C) (WMLCIA), I A=MXB1 ,MXBC) 

RETURN 

1C2C FORMAT ( IX, 8G 16. 7 1 

1530 FORMAT l 2FP T , 5GX , 16H STRE AMLl NE PLOTS ) 

1 54C FORMAT < 1 X , 7G 1 6 . 7 / ( 1 SX , 6G 1 8 . 7) ) 

155C FORMAT ( IF. 1 , 26X , 22HSTREA ML I NE COCRDI NATES///7 X ,8HM CCCKC., 
1 3( SX, 1CHS1R EAM FN.,SX,bH THETA )//) 
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156C FORMAT { 2HPL , 4GX , 7 CH STREA ML! NE S ARE FLCTTEO WITH ThETA ACRCSS THE 
1PAGE AND M DOWN THE PAGE) 

160C FORMAT (1FC.22HWM ARRAY (RHO* W- SUB-M ) ) 

16IC FORMAT 1 1FG, 40HWMU ARRAY (RHG*W-SUB-M CN UPPER SURFACE)) 

162C FORMAT I 1FC, 4ChWML ARRAY (RHO*W-SUa-M CN LLwER SURFACE)) 

END 


SUBROUTINE TASVEL 

REAL K ,K1,K2,K2,K4,K5,LMAX,LMI N , LAMBDA , MU , ML , M R , MXU ,MX L » MPL 
INTEGER SRW, FIR ST .CASE ,6LDA TA ,E RPR T , S T RFN , S LCRD ,SLPLT > ARPRT ,SURVEL 
COMMON SRW. ITER , GAM , AR , T I P , RHOI P , WTFL ,CMEG A ,L AMBDA , CP , EXPON , PIT CH, 

1 CHORD, STGR.BETAI ,BE TAO.D TLR ,RI ,ALUI ,ALLI ,RC , A LUC , ALLC , 

2 MXB I,MXBO,MX,NBBI , NUSP , NL SP , NR SP , M NT , V TOL , 

2 BLCATA,NULAKI,ERPRT,STRFN ,SLCRD , ARPRT, I NTVEL .SURVEL, 

4 MOI 50) ,XSPU( 50) ,ML (50) , XSPL (50) , MR 150) , RMS P 150) , BES P 150 ) , 

5 W, WR, TOLER, BDA , BOO, U 125001 , A 12500, 4) ,X 12500) , PHD 12500), 

6 DXC2U1 100 J ,OXDZL 1 100) .SLOPE (50 , EMU 150) , SLL PE 150 ) .EML150), 

7 RMl ICO), BE! ICC) , SAL 1 IOC ) , XL 1 IOC ) , XL (100 ) , RMU 1100 ) , RML 1 100 ) , 
a NUl 100 ) ,NL ( ICC) ,UINT( 11) , XI NT 111) , MPL 1100) , 

S HA, FB, NXN,MXBIMl,MXBOPl,JU,JL,HU,HL,NBdC,NBUC,NCH, 

1 18TE1, IBTE2, ITE2.HA1TE .HANTE ,UNI , UNC ,T W W , I T E P A , 

2 KHOLl ICC) ,RHOL( 100) .RHOUT11CO) .RHCLTUOO) .BEUUOO ) ,B ELI 100), 

3 AAA (IOC) 

C THE FULL CW ING VARIABLES ARE ALL IN CCMMCN BY EQUIVALENCE 

DIMENSION WM ( 2 500 ) , WX 1 2500 ) ,V(2500) , BETA 12500 ) ,SL 1 1 100 ) 

DIMENS ION wMUl ICO) , WML 1 100) ,WXU (100) ,W XL (100) ,MXU( 100 ) ,MXL(100 ), 

1 BETAU(IOO), BETALl ICO) ,WL(10C) ,WL (100) , 

2 USP(ICC), IU1 lCO.ILUUC) 

EQUIVALENCE (A(1,1),WM(1)),(A(1,2),WX(1)),IA(1,3),V(1)), 

1 (All, A), BETA! 1) ), (K( 1401) ,SL(1) ) 

EQUIVALENCE l K ( 1 ) ,WMU( 1) ),(K(1C1),WML(1)),(K(201),WXU(1)), 

1 (K ( 2C1 ) ,WXL ( 1 ) ) , ( K( 401 ) ,MXU (1) ) , 1 K (501 ) ,MXL(1> ) , 

2 (K1CC1), BE TAU ( 1) ) , ( K ( 70 1 ) , BE TAL 1 1 ) ) , 1 K (80 1 ) , WU ( 1 ) ) , 

3 (K( SCI >,WL( 1)),(K( 1001) ,0SP(1)) ,1X11101) ,IU(1)) , 

A (K( 1201), 1L(1)) 

UIMENS ION QTDMU ICC) ,D TDMLUOO) ,BBB (100) 

DIMENS JON XDO WN 14CC),YACR0S(40C) ,KKK(14) 

EQUIVALENCE 1 K 1 130 1 ) ,D TO MU ( 1 ) ) , l K (1 401 ) ,0 TDML 1 1 ) ) , 

1 (K( 1501), XDOWNl 1) ) , 1X11 SCI) .YACRCSU) ) 

C 

C CALCULATE R HO *W— SUB- THETA 
C 

C START AT 1A = 1 (CASES 1,2,3) 

C 

CASE =J 
JU = C 
JL = C 
KK 1 = C 
KN = C 
IB = C 
I A 1 = 1 
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hAl = HA 
IAN = MXB I 
I = N B E I + I 

ENC UN UPPER SURFACE (CASE 1) 

£1C IF( IB .CT.NLC 1) > GU TO 600 
DU 53 C IA = IAN » M XBU 
53C IF(NU( 1AI.GT.IB) GU TO 540 
CASE = 2 
GU TO £50 
54C IAN = 1A 

I A = I AN- I 

XI = FLQA T( IB J*hB 

CALL 1NTPL 1( XU, IA , XI ,HAN , 1 , K4 , MU , XSP l , ALUI , A LUG , NUS P,S LU PE , EMU , 

I HI , RO, CHORD, STGR ,PI TCH ,1. ,HA ,HB , PXBI ,MXBO , RM { MXBI ) , RM( MXBU) ) 
GO TU 770 

ENC AT IA = MX (CASE 2) 

55C I F ( IB.NE. IBTEI ) GO TO 554 
IAN = MXBU 
hAN = PANTE 
KN = 1 
GO TO 7 7 C 

554 CO 555 IA=MXBI,MXBO 

555 IF( IB .C-T.NLI IA ) ) GO TO 560 
IF( IB .GT.NL( 1)) GO TU 6GC 
IAN = MX 

hAN = PA 
GO TO 77C 

ENC ON LONER SURFACE (CASE 31 
56C CASE = 3 

I F( IB .GT.NL( 1 ) ) GU TO 6CC 
DO 57C I A =MXB I , M XBQ 
57 C IF(NL( IA) .LT.IB I GO TO 580 
I A=MX BO 
58C IAN = IA 
IA = IAN-1 
XI = FLOAT! IB )*HB 

CALL INTPL 1( XL , IA , XI, HAN, I , K4 , M L , XS P L , ALL I , AL LO , NLS P ,S L LPE , EML , 

1 R I , RO, CPORO, STGR ,P I TCH ,-l. ,HA,HB ,MXBI , MXBO , RM ( MXB I) , RM ( MX BO ) ) 
KN = 1 
GU TO 77C 

START ON LOWER SURFACE (CASE 41 
60C CASE = 4 
KK 1 = 1 

62C CO 63 C IA =MXB I ,M XBO 
6 3C IF(NL( 1A+ 1 ) . GT .NL ( IA ) ) GO TO 64C 
CAS E= 5 

go to teo 

6 4C IB = NL( IA )+l 
I A 1 = IA 

65C DO 66 C IA = IA I , M XBO 
66 C IF(NL( 1A+17.GE.IB I GO TO 67C 
CASE = 5 
GO TU ECO 
670 IA1 = IA 
IA = IA1+1 
XI = FLOAT (IB )*HB 
I = IL ( IA ) + IB-NL( IA ) 



CALL IMPLH XL , I A , X 1 f HA 1 , 0 , K3 , H L , XSP L , ALLI ,A LLC , NLS P , S LLPE, EML, 

1 K I , RO, CHORD » STGR *PI TCH,— 1. ,HA,HB,MXBI , MXBO , RM ( MX 81 ) , KM (MX BO) ) 
GO TO 74C 

C START ON UPPER SURFACE (CASE 5) 

ESC DO 69C I A =MXB I »M XBO 
69C IFCNUC 1A+1 I.LT.NUUA) ) GO TO 7CC 
CASE = 6 
GO TO 765 

70C IB = NUC I A )— 1 

I A I = IA 
KK 1 = 0 

71C CU 72 C IA = IA 1 , M XBO 
72 C 1FCNUC IA+1 I.LE.IB ) GO TO 73C 
CASE = 6 
GO TO 765 
73C I A 1 = IA 
IA = IA +1 
XI = FLOA TC I B )*HB 
IFC IB .EQ.NUCMX )) X1=STGR 
I = I U ( IA HIB-NLC IA ) 

CALL INTPL 1C XU ,IA ,X1 ,HA1 , 0 , K3 , PC , X SPU , ALUI ,ALUC,NUSP,S LUPE.EMU, 

1 k I, RO, CHORD, STGR, PI TCH, 1. ,HA,HB,MXBI ,MXBQ » RK ( MXB I ) , RMC MX BU ) ) 

C CASES 4 AND 5 
74C I AN = IA I + 1 

00 75 C IA = IA N, MXBO 

75C IFCNLC IA).LT.IB) GO TO 76C 
IFC IB.NE. IBTE1 ) GU TO 755 
IAN = MXBO 
HAN = FANTE 
KN = 1 
GO TO 7 70 
C ENC AT IA = MX 
755 IAN = MX 
hAN = FA 
KN = C 
Go TO 77C 

C ENC ON LOKER SURFACE 
76C IAN = IA 
IA = l A— 1 
KN = 1 

CALL IN TP Ll( XL , IA ,X1 ,HAN, 1 ,K4,ML,XSPL,ALLI ,ALLC,NLSP, SLOPE, EML , 

1 R I , RO , CHORD, STGR, PI TCH ,-l. ,HA ,H B , MX BI , MXBO , RM C MXBI ), RMC MXBO) ) 
GO TO 77C 

76 5 IFC IBTE2.EQ. 1CCC) GO TO 78C 
I A 1 = MXBO— 1 
HA 1 = FA1TE 

1 = ITE2 
IAN = MX 
HAN = FA 
IB = IETE2 
KK 1 =1 

77C CALL VELOC ( U , MPL , IA 1 , I AN ,H A1 ,H AN , I ,1 B , MX , NXN ,NBBO , JU , J L , KK 1 , KN, 

I USP ,WX,WXU,MXU,WXL,MXL,NU,NL,AAA> 

18 = ie+i 

I = 1+) 

IFC CASE. EQ. 5) IB = IB-2 

GO TO (51C,55L,56C, 650, 710, 780) .CASE 
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78C 


CALL SORT X Y (MXL,WXL,JL) 

CALL SPLINT (MR ,8E SP ,NR SP,KXL,JL,8EL»AAA) 


CALL SPLINT (MR ,BESP ,NRSP ,MXL ,JL,BEL,AAA) 
CAlL SPLINT (MR ,RMSP, NRSP, MXU,jL,RMU, AAA) 
CALL SPLINT (MR fRMSP ,NRSP »MXL,JL»RML»AAA) 
LAST = C 
CU 79C IA = 1 , M X 
FIRST = LAST+1 


LAST = F IR ST +NL ( lA)-NU(IA) 

00 79C I = F IR ST »LA ST 
79C WX( I ) = -WXI I)/BE(IA J*wTFL 
00 800 1=1, JL 

£00 WXU( I ) = — WXU ( I)/BEL(I )* WTF L 
DO 81C 1= I, JL 

81C MXL ( I ) = -WXL( I )/dEL< I )*WTFL 


ENC OF R HG*W- SUB— T FETA CALCULATION 

IF( ARPRT.LE.C ) GO TO £3C 
WRITE (6,1630 
L AST = 0 

DO £2 C I A = 1 , M X 
F IkST = LAST+ 1 

LAST = F 1RST+NL! IA )-NU( IA) 

82C WRITE! 6, 1020 ) (MX(I), I =F I R ST , LA S T > 

WRITE (6,‘l64C) 

WRITE! 6, 1C 20) (MXC( I ) , WXU (1 ) , 1=1, JL) 

WR ITE (6, 1650 

WRITE( 6, 1C20 ) (MXL( I ) , W XL ( I ) ,1=1, JL) 

CALCULATE RFO*W AND ANGLES AT INTERIOR POINTS 

€3C CON T IN LE 
LAST = C 
CO 850 IA = 1, MX 
FIRST = LAST+1 
LAST = FIRST+NLI IA )-NU(IA) 

DO 850 I = FIKST,LAST 

V(I) = SWRT( WX( I )**2+WM( I )**2) 

IF(WM( D.EQ.C. ) GO TO £ AG 
BETAl I ) = A TAN (WX(I ) / WM ( I ) )*57. 29577 

GO TO £50 
£AC 8 E TA( I 1 = 90. 

£50 CONTINCE 

IF( ARPRT.LE.C) GO TO 87C 
WRITE (6, 1660 
L AST = C 

CO 86 C IA = 1 , M X 
F IK ST = LAS T + 1 

LAST = F 1RST+NL( lA)-NL(IA) 

860 WRITE! 6, 1C2C ) ( V ( I ) , I =F I R S T ,LA S T) 

CALCULATE OENSITY AND VELOCITY AT EACH PCINT 


87C LAST = C 
R EL ER = 0 . 

IF( IN TV EL ,GT .0 ) WRITE (6,1730 
CP T IP = 2 .*CP* T IP 
CO 5 1 C IA = 1, MX 
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FIRST = LAST+1 

LAST = FIRST+NLdA )-NU(lA) 

TwLMR = 2 .*OMEGA*LAMbDA-(UMEGA*RM(I A) ) **Z 
DU 9GC I = F IR ST »LAST 
RHUT = RHO { I J 

CALL CENSTY (VU) ,RHU ( 1 1 , VEL , TWLMR ,CP T I P * E XPON , RHO I P ,G AM, AR » T IP* 

1 VTUL ) 

RELER = AMAX1(RELER,ABS( lRHCT-RHO(I ) ) /RHC(I) ) ) 

SOC V( l ) = VEL 

IF! IwTVEL .LE.C) GU TU 91 C 

WRITE! 6, 1 6 7 C ) I A » (VU), faETA(l) ,1 =FIRST » LAST ) 

SIC CONTINUE 

IF! ITER.GT.G ) ITERA = ITER 
WRITE (6,1735) ITERA, RELER 
ITtRA = ITER+1 

IF! RELER.LT. .CC1. AND. I TER. GT.l ) ITER = -1 

IF! GAM ,EQ .1. 5 .ANO.AR.Ey. 1GOO. .AND. TI P. EC. 1 . £6 ) I TER=G 
CALCULATE SURFACE VELOCITIES BASED CN AXIAL COMPONENTS 

IF! SURVEL .GT .C ) WRITE (6,1660 
BETAU! MXB I ) = SC. 

BETAL ( MXB I ) = -SO. 

WU(MXB 1 ) = 0 . 

WLIMXB 1 ) = 0 . 

BETAU! MXBO ) = -SC. 

BETAL(MXBO) = SC. 

WU(MXBO) = 0. 

wl(mxbc) = o. 

MX B IP 1 =MXB 1+1 
MX BOM 1 =MXBQ- 1 
AAA(MXEI) = C. 

BBE(MXEI) = C. 

DO S2C 1A=MXBIP1,MXBUM1 

AAA! I A )=AAA( IA- 1 I+SURT (HA**2+((XU(IA)~X01 I A— 1 ))*(RM(IA)+RM(IA-1)) 

1 / 2 . )** 2 ) 

TANTHU=CXC2U( IA 1*RM(IA ) 

BETAU! IA) =A TAN! TAN THU )# 57. 2S577 
WMU! IA )=WMU( IA )*SWRT( 1 . + TA N IH 0* TANTH C ) 

TWLMR = 2 .*UMEGA*LAM6DA-(0MEGA*RM( I A) ) **Z 

CALL CENST Y (WML! I A ) ,RHOO( I A ) , WOd A) , T WLMR ,C PT I F , EX PON , RHC I P , GAM, 

1 AR , T IP , VTOL ) 

BBB! IA )=BBB( IA- 1 ) + SQR T (HA* *2+ ( (XL (I A) -XL! IA-1 ) ) * ( RM ( I A ) +RM! IA-1) ) 

1 /2.)**2) 

TANThL =DXDZL ( I A )*RM( I A ) 

BETAL! IA) -A TAN! TANTHL)*5S. 2S577 
WML! IA ) = WML ( IA )* SQR T ( 1 . + TA N THL* TANTHL) 

CALL CENSTY ( WML ! IA ) ,RHOL ! I A ) , W L !I A ) , TWLMR ,C PTI F ,EXPON, RHCI P, GAM, 

1 AR , T IP , VTUL ) 

S2C CONT IN LE 

IF! SURVEL .LE .C ) GO TO S27 

A A A ! M X BO ) =AA A ( M XBOM 1 ) + SORT (H A* *2+ ( (XL (MXBO) -XU! MXBCM1) ) *( RM( MX 80) 

1 +RM(MXBUM1) )/2.)**2) 

BBB(MXEO)=rtBB(MXBOMl)+SORT(HAv<‘2+!(>L(MXBC)-XL(MXBGMl)) : *(RM(MXBO) 

1 +RM ( MXBOM l))/2.)**2) 

WRITE (6, 1 6SC ) (MPL! IA) ,WC(I A) , BE T A L ( I A ) , AAA ( I A ) , WMU ( I A ) , W L ( I A ) , 
1 BETAL! IA ),BBB(IA) , wML II A ) ,1 A=MXBl ,MXBO) 

NP1 =C 

DO S23 IA =M XB I ,M XBQ 
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IF1ABSIBETAUI IA ) J.GT.6C.) GC TG 923 
NP1 = NP141 

Y ACROS (NP 1 ) = kiUIA) 

XDUWN I NP 1 ) = MPL ( IA) 

92 2 CONTINUE 
NP2 = NP1 

DO 925 IA=MXBI,MXBO 

IFI ABSI8ETALI IA ) ) .GT.60. ) GO TO 925 
NP2 = NP2-U 

Y ACRUS (NP 2 ) = WLIIA) 

X DOWN ( NP 2 ) = MPL(IA) 

92 5 CONTINLE 
NP3 = NP 2 
NP2 = NP 2— NP 1 

CALCULATE SLRFACE VELOCITIES BASED ON TANGENTIAL CCMPGNENTS 

WRITE (6, 1700 
927 CONTINLE 

CALL BLDDE 1(MU»XSPL» SLUPE ,EMU,NUSP,RI ,ALUI »RG , ALUO »CHORC»ST GR» 

I P ITCH, 1., JU,MXU,MXBI ,HA , D TO PU , RM ( MXBI ) ,RMt MXBQ) ) 

CALL BLDDE 1(HL ,XSPL, SLLPE .EML.NLSP ,RI , ALLI , RO , ALLO , CHGR C,ST GR, 

1 PITCH,— 1 . , J L , M XL ,HXBI ,HA , D TDML ,RM ( P XB I ) ,RM1MXBG) ) 

C UPPER SURFACE 

BETAUI 1) = 9C . 

WXU( I ) =AB S ( WXU ( 1) ) 

TWLMR = 2 .*0M EGA* LAMB DA— IUMEGA*RMU ( 1 ) ) **2 

CALL CENSTYIWXUt 1),RH0UT( 1),WU I 1 ) ,TWL HR , C PT I P , EXPCN , RHOI P, 
1 GAM ,AR ,TIP, VIOL) 

BETAUI JU ) = — SC • 

WXU(JU)=ABS(WXUIJUn 

TWLMK= 2 .*QM EGA* LAMBDA— IQMEGA*RMUIJU) ) **2 

CALL CEN STY ( WXU I JU ) , RHOUT ( JU) »WU I JU) , TWL PR ,C FTI P , EXPON, RHO I P, 
1 GAM.AR.TIP, VTOL) 

J UM 1= J L— 1 

DO 540 I = 2 , J UM I 

TANTHU = OTOMUl I)*RMUII J 

BETAUI I) =A TAN l TAN THU)*57. 25577 

IFI TANTHU.EQ.O.) GO TO 530 

WXUI I )= ABSt WXU I I ) )* SQR TI l. + l. /ITAN1HL*TANTHU) ) 

TWLMR = 2 .*OMEGA* LAMBDA— I OMEGA* R MU 1 1 ) ) **2 

CALL CENSTY I WXU 1 1 ) ,RHOU T 1 1 ) ,WLII) ,T WLMR ,C PTI P , EXPGN, RHC I P , GAM, 
I AR, TIP, VTOL) 

GO TO 540 
53C WUl I) =0. 

54 C CONTINUE 

IFISURVEL .LE .0 ) GO TO 547 

WRITE 16,1710 IMXLI I ) , WU 1 1 ) , BETAU II ) ,WXU II ) , 1=1, JU ) 

00 545 I = 1 , J U 

IFI ABSIBETAUI I ) ).LT.3G.) GO TO 545 
NP3 = NP3+1 
YACR0SINP3) = WUII) 

X00WNINP3) = MXLI I ) 

545 CONTINUE 
NP 4 = NP 3 
NP 3 = NP 3— NP 2— NP 1 
C LOWER SURFACE 

WRITE 16,1720 
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94? CONT IN LE 

CO 96C 1=1, JL 

TANTHL = D TOML ( 1)*RML(I) 

BETAL(l) =A TAN ( TANTHL 1*57. 29577 
IF { TANTHL .EQ.C . ) GO TO 950 

WXL ( I ) =AB Si wXL ( 111* SQR T( l.+l./!TANTHL*TANTHL) 1 
TWLMR = 2 .*UMEGA* LAMBDA— (OMEGA *RML (111 **2 

CALL CENSTY ( wXL ( I 1 ,RHOLT U 1 , WL ( 1 1 , TWLMR ,C FT 1 F , EXPON, RHGI P, GAM, 

1 AR , T IP , VTOL 1 
GO TO 96C 
S5C WL( II =0. 

56C CONTINUE 

IF ( SURVEL .LE .0 ) RETURN y 

WRITE (6.171C1 IMXL( I 1 , WL(I 1 , BETA L( 1 1 ,W XL U 1 , 1=1, JL) 

00 970 1=1, JL 

IF( ABS(B£TAL( i 1 1.LT.30.) GO TO 970 

NP4 = NP4+1 

YACRO S (NP 4 1 = WL ( I ) 

XDOWN( NP4 ) = MXL { I 1 
97 C CONTINUE 

NP4 = NP4-NP3-NP2-NP1 

KKK( 11=0 

KKKI21 = 4 

KKR ( 3 1 = NP1 

KKKC 5 1 = NP 2 

KKK, { 7 1 = NP 3 

KKM 9 1 = NP 4 

P = 5 • 

WRITE (6, 1740 

CALL PLGTM Y( XDGWN , YACROS *KKK ,P) 

WRITE (6,1750 
RETURN 

1020 FORMAT (1X.8G16.7) 

163C FORMAT (1F1,26HWX ARRA Y (RHO* W- SUB— THE I A) 1 

1640 FORMAT ( IhK , 4( 6X , 3HMXU , 13X , 3H WXU ,7X1 /48H <M CCLRC. VS. RHO*W-SUB-T 
1HETA ON UPPER SURFACE) 1 

165C FORMAT ( 1 FK, 4 ( 6 X , 3HM XL , 13X , 3H WXL ,7X) /48H (M CCCKC. VS. RHO*W-SUB-T 
1HETA UN LOWER SURFACE) ) 

166C FORMAT ( 34F1ARRAY OF RHO*W A1 INTERIOR PCINTS ) 

167C FORMAT ( 1FL , 3H IA = , I 2, 5 ( 24H VELOCITY ANG LE ( DEG) )/ ( 3X , 

1 5<G15.4,F9.2))1 

168C FORMAT ( 1K1 , 15X , 1H* ,2 IX , 44H SURF ACE VELOCITIES BASED ON AXIAL COMPO 
IN ENTS , 45X , 1H* , / 16 X, 1H* , 12 X , 16H UPPER SURFACE »25X , 1H*, 15X , 16HL0WE 

2R SLR FACE, 25X.1H* ,/7X,IHM,8X ,1H* ,2 (3X,8HVELCCIT Y ,3X ,43H ANGLE! DE 

3G) SURF. LENGTH RHG*W * )) 

1690 FORMAT (IF ,G13.4,3H * ,G 12 . 4 ,F 9. 2 ,2G1 5.4 ,faH * ,G12.4,F9.2, 

1 2G1 5.4, 3H * ) 

170C FORMAT (51H1 SURFACE VELOCITIES BASED CN TANGENTIAL COMPONENTS/ 

1 25X , 13HUP PER SURF AC E/7X ,1HM,10X,8HVELCC1TY ,3X ,10 H ANGLE! OEG) , 

2 3X , 5HR FO* W ) 

171C FORMAT (IF , 2G 1 3 . 4 ,F 9. 2 ,G 1 5. 4) 

172C FORMAT ( 2 5X , 13FL0WER SURFACE/7X ,1HM ,1 OX ,8H VELCCITY ,3X , 

1 lOFANGLE(DEG), 3X,5HRHO*W) 

1 73 C FORMAT ( IF 1, 4CX, 34HVELOCI TIE S AT INTERIOR MESH PCINTS ) 

1735 FORMAT ( 14HL ITERATION NO., I 3 , 3 X ,36H MAXI MUM RELATIVE CHANGE IN DEN 
1SITY = ,G1 1 .4 ) 
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174C FORMAT < 2FP T , 5CX , 24H8 LADE SURFACE VELOCITIES) 

175C FORMAT ( 2HPL ,3 7 X , 82H VE LGC I T Y (ME TER S /SECCND) VS. MERIDIONAL STREAM 
1LINE 0 ISTANCE(METERS) DOWN THE PAGE /2HPL/ 

2 2HPL.5CX, 43H* - UPPER SURFACE, BASED CN AXIAL COMPONENT / 

3 2HPL.5CX, 43H+ - LOWER SURFACE, BASED ON AXIAL COMPONENT / 

A 2HPL.50X, 48HC - LPPER SURFACE, BASED CN TANGENTIAL COMPONENT/ 

5 2HPL.5CX, A8HX - LOWER SURFACE, BASED CN TANGENTIAL COMPONENT) 

END 


Subroutine DENSTY 


■' DENSTY calculates the subsonic relative velocity W and corresponding density p 

that result in a given value of the mass flow parameter pW. This is done using equations 
(B5) and (B6) which are an algorithm based on Newton’s method, 
j If the value of pW is too large, there is no solution. In this case an error message 

I is printed out. If a ’’continue" control card is used (see p. 86), W cr and the corre- 

j sponding density are calculated as output and the program continues. This makes it pos- 

| sible to get an approximate solution even though there may be one or two points with too 

large a value for pW. 

I The input arguments are as follows: 

f RHOW pW 

! 

i RHO initial estimate for p(p- n may be used) 

i TWLMR 2coA - (wr) 2 

[ CPTIP 2cT! 

\ P in 

j EXPON l/(y - 1) 

u 

RHOIP p! 

j m 

I GAM y 

f AR R 

TIP T? n 

i 

I VTOL convergence tolerance on relative change in W 

I 

S' 

| The output arguments are as follows: 

I RHO p 

r. 

I VEL W 

I 

j The internal variables are as follows: 

; RHOT newly calculated estimate for p 

| RHOWP d(pW)/dp 


i 
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TEMP 


(2-y)/(y-l) 


(T/T! n ) 

TGROG 2yR/(y + 1) 

TTIP t / t : 

in 

VELNEW newly calculated estimate for W 


SUBROUTINE DENSTYIRHOW.RHO » VEL , TWLMR ,C PTI P ,EX PCN ,RHCI P , GAM, AR, T IP, 
1 VTOL) 

VEL = RFOW/RHO 

10 TTIP = 1 I VEL** 2+TWLMR ) /C PTI P 
1 F( TT IP.LT.O. ) GU TO 30 
TEMP = TTIP**(EXPON-l.) 

RHOT = RHO IP*TEMP*TTIP 

RHOWP = -VEL**2/GAM*RH0IP/AR*TEMP/TIP«-RH0T 
IFtRHOUP.LE.C. ) GO TO 30 
VELNEW = VEL + ( R HOW- RHOT* VEL) / RHOWP 
I FI ABS( VELNEW- VEL l/VELNE W. LT- VTOL) GO TO 20 
VEL = VELNEW 
GO TO 10 
20 VEL = VELNEW 
RHO = RHOW/VEL 
RETURN 

30 CALL ARERK < 29H VALUE OF RHO*W IS TOO LARGE*) 

TGROG = 2.*GAM*AR/(GAM+1.) 

VEL = SQRT ( TGRUG*TIP*(1.-TWLMR/CPTI P) ) 

RHO = RHO IP* < 1 { VEL**2+TWLMR) /CPTI P) **EXPGN 

RETURN 

ENC 


Subroutine BLDCR1 

BLDCRl obtains the 6 -coordinates and the slopes of the blade surfaces correspond- 
ing to the given m-coordinates. BLDCRl may be used in two ways, either to obtain the 
information at all vertical mesh lines in one call, as when called by COEF directly, or at 
a single specified point, as when called by INTPL1. The value of NCH is the number of 
points at which output is desired. Since BLDCRl is used for either the upper or lower 
blade surface, SURF is used as a code to determine which surface is desired. SURF = 1. 
for the upper surface and SURF = -1. for the lower surface. 

The entire blade surface is defined by the leading- and trailing- edge radii and by two 
cubic spline curves (upper and lower surfaces), which are piecewise cubic polynomials. 
The procedure then is to scan the spline points to determine which interval the 
m-coordinate (ZINT) is in, and then to calculate the 0-coordinate and derivative, both of 
which are specified analytically. 
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The input arguments are as follows: 

Z array of m-coordinates of spline points for blade surface (upper or lower) 

XSP array of 9 -coordinate of spline points for blade surface (upper or lower) 

SLOPE array of slopes at spline points for blade surface (upper or lower) 

EM array of second derivatives at spline points for blade surface (upper or lower) 

NSP number of spline points on blade surface (upper or lower) 

RI see fig. 10 (p. 15) 

ALI either ALUI or ALLI (see fig. 10) 

RO see fig. 10 

ALO either ALUO or ALLO (see fig. 10) 

CHORD see fig. 10 

STGR see fig. 10 

PITCH see fig. 10 

SURF code to indicate upper or lower surface, SURF = 1. , for upper surface and, 
SURF = -1., for lower surface 

NCH number of points for which output is desired 

ZDSTT used as input only when NCH = 1, then it is m -coordinate for which correspond- 
ing 9 -coordinate for blade surface is desired 

MXBI same as main program, number of mesh points on line AB 

RMI r at leading edge 

RMO r at trailing edge 

The output arguments are as follows: 

X array of 9 -coordinates at the vertical mesh lines for blade surface, or if 

NCH = 1, at m = ZINT 

DXDZ array of slopes at same points as X 
The internal variables are as follows: 

HA basic mesh spacing in meridional (m) direction 

LA index of vertical mesh line 

IFST index of first point considered 
ILST index of last point considered 
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K index of spline point 

RMZ difference between m-coordinate of point considered and m-coordinate of center 
of leading- or trailing-edge radii 

SRW integer variable in common used to obtain output useful in debugging; when 
SRW = 19, BLDCR1 will write out calculated blade coordinates and corre- 
sponding slopes 

SW coefficient with value of zero on upper blade surface and 1 on lower blade sur- 

face; used to add pitch to computed blade coordinate for lower surface only 

ZINT m-coordinate at which Q -coordinate and slope of blade surface are required 


SU6R0UTIN 6 BLOCR 1 1 Z , XSP. SLOPE ,E M.NSP ,RI »ALI , RC > ALC , CHORD , S T GR, 

1 PITCH, SUR F , X ,NCH ,ZINT ,MXBI ,DXQZ, RMI ,RMO) 

C 

C SURF =1. — UPPER SURFACE 

C SURF = -1. — LOWER SURFACE 

C 

OIMENS ION XSPINSP ) , Z ( N SP ) ,X(NCH) ,SLCPE(NSP) ,E Ml NSP) ,DXDZ I NCH) 
COMMON Skw 
INTEGER SRW 
SW = G. 

IF( SURF.LT.O.) SW = 1. 

IFST = 1 
ILST = I 

IFINCF.EQ.I) GO TU 10 

IFST = MX B I 

ILST = NCF+MXB I- 1 

HA = CFORD /FLOAT! NCH- 1 ) 

ZINT = C. 

1C K = 2 

00 IOC IA = IF ST , ILST 
20 IFIZ INT.GT.ZI 1)) GO TO 30 

XI IA) = SQR T ( Z IN T* I 2.*RI-ZINT) ) /RMI *SURF+PITCH*SW 
RMZ = K I- ZINT 

IF! IA.NE. IFST) OXOZ(IA) = RMZ/SQRT ( RI **2 — RMZ**2 ) *SURF/ RMI 
Z INT = ZINT + HA 
GO TO ICO 

30 IF(ZINT.LE.ZIK) ) GU TO 50 
IF(K.GE.NSP) GO TO 60 
K = K + I 
GU TO 2C 

5C S = Z I K)-Z(K- i ) 

X( IA ) = tM(K-l)*(Z<K)-ZINT)**3/6./S+EMlK)*(ZINT-Z(K-l) )**3/6./S 

1 + ( XSPI K )/S-EM(K)*S/6. )* I ZI NT- Z IK-1 ) ) + (XSPIK-i ) / S -EM ( K-l ) *S/ 6.) 

2 *( Z (K )-Z IN T ) 

OXOZ ( I A ) = -EM I K-l)* < ZIKJ-ZI NT) **2 /2 . /S+EMIK) *(Z(K-1)-Z INT ) **2/2. 
1 /S+IXSPIK >- XSP (K-l) )/S- I EM (K) -EM l K-l) > *S/6. 

ZINT = ZINT+HA 
GU TO IOC 

6C IF I I IA.EQ.NCH) .AND .(1A.GT.1) ) GO TO 70 

X( IA) = STGR + SURF* SQR T (1C HURU-ZI NT) * (2 . *RO-CHORD+Z I NT ) )/RM0 + PITCH 
I *Sw 

RMZ = CHORD-ZINT-RQ 
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IF( IA.NE. ILST) DXOZ(IA) = RMZ/SQRT ( RC**2— RMZ**2 ) *SURF/ RMO 
ZINT = ZINT+HA 
GO TO ICG 

70 X( IA) = STGR+P ITCH*SW 
100 CONTINUE 

IF( SKW .EO . 19 ) WR I TE ( 6 » 1GG0 ) { X ( I A) , 0 XDZ ( I A) , 1 A= I FST , I LST ) 

RETURN 

1000 FORMAT ( 1 X , 5 4H IN TERPOLA TED COORDINATES AND SLOPES CCMPUTEC EY BLUC 
1RD, /< 5X f 2 E 16 .8 ) ) 

END 


Subroutine BLDDE1 

BLDDEl obtains the slopes of the blade at given m -coordinates. It is used by 
TASVEL to obtain the blade slopes at each horizontal mesh line. BLDDEl is similar to 
BLDCR1, except that the m -coordinates are an input array and the 9 -coordinates are not 
given as output. 

The input arguments for BLDDEl are the same as those for BCDCR1, except that 
ZINT is not input. Also included are 

ZX array of m -coordinates from the line BG in fig. 4 for which slopes for blade 

surface are desired; these values must be arranged in increasing order 

HA basic mesh spacing in axial direction 

The output of BLDDEl is 

DXDZ array of slopes at m -coordinates in array ZX 

The internal variables are as follows : 

IA index of point in ZX array 

K index of spline point 

RMZ difference between m-coordinate of point considered and m-coordinate of center 
of leading- or trailing-edge radii 

SRW integer variable in COMMON used to obtain output useful in debugging; when 
SRW = 20, BLDDEl writes out calculated blade slopes 

ZINT m-coordinate from blade leading edge at which blade slope is desired 
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SUBROUTINE BLODE 1 tZ ,XSP , SLUffc ,EM,NSP,RI,ALI , RC , ALC , CHORC,ST GR, 

1 P ITCH, SURF ,NCH,ZX,MXBI ,HA,DXUZ ,RM ,RMG) 

C 

C SURF =1. — UPPER SURFACE 

C SURF = -J. — LOWER SURFACE 

C 

DIMENSION XSP(NSP)»Z(NSP) .SLCPEINSP) ,E M (NSP) , DXQZ ( NCH ) » ZX ( NCH ) 
COMMON SR W 
INTEGER SR W 
1C K = 2 

CO IOC IA = 1»NCH 
Z INT = ZX( IA ) 

2C IFIZlNl.GT.Ziin GO TO 30 
RMZ = RI-ZINT 

IF< IA.NE. l.OR.SURF.LT.C.) OXOZIIA) ■= RMZ/SCRT ( R I -RMZ **2 ) *SU RF 
I /RM I 

GO TO ICG 

30 IF(ZINT.LE.Z(K ) ) GO TO 50 
IF(K.GE.NSP) GO TO 60 
K = K + I 
GO TO 20 

5C S - Z ( K )— Z ( K— 1 ) 

DXDZ( I A ) = -EMIK-l )* ( Z(K)-ZI NT) **2 /2. /S+EMIK) * (Z IK-1 )-ZINT )**2/2. 

1 /S+IXSPIK )— XSP ( K— 1 ) ) /S- ( EM K) -ERIK-1)) *$/6. 

GO TO ICO 

6C RMZ = CFORO-ZINT-RO 

IF I ( I A. NE.l. AND. IA.NE. NCH). UR. SURF. LT.O.) DXOZ ( IA) = RMZ/ SORT ( R0**2 
1-RMZ ** 2 )* SUR F /RMU 
GO TO ICO 
IOC CONTINUE 

IF( SRW .EQ .20) WRI TE 16, 1C00) (ZX(IA) ,DXDZ(IA),IA-1,NCH) 

RETURN 

1C0C FORMAT (lX,5fchZ CCORO. AND INTERPOLATED DERIVATIVES COMPUTED BY BL 
1CDER, /(5X, 2E16.B)) 

ENC 


Subroutine INTPL1 

To compute the terms of the matrix A of equation (A7), it is necessary to obtain the 
distance along a horizontal mesh line from a mesh point near the blade to the blade itself. 
This is the quantity h 3 or h 4 in equation (Al). This value is computed by INTPL1. 
Since the equation of the blade surface is known, this amounts to finding the root of an 
equation. The root is found by INTPL1 by an iterative procedure, sometimes called the 
method of false position (falsi reguli) . The variables shown in figure 16 correspond to 
those used in INTPL1. After H has been calculated, the actual value of the spline curve 
(XI) is computed by BLDCR1 and a reduced interval is considered so that the curve still 
crosses the value XI. Then the procedure is repeated on the smaller interval. A few 
iterations will determine the value of z for which the spline curve crosses the mesh 
line, and from this, h 3 or h 4 is determined. 
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The input arguments are as follows: 

HA basic mesh spacing in axial direction 

HB basic mesh spacing in blade-to-blade direction 

IA index of vertical mesh line on which mesh point lies 

MXBI number of mesh points on line AB 

N integer which is zero when hg is to be computed and which is 1 when h 4 is 

to be computed 

X array of blade 0 -coordinates at each vertical mesh line 

XI 0 -coordinate of mesh point considered 

The remaining input arguments are transmitted to BLDCR1. Their definitions are in- 
cluded in the description of this subroutine. 

The output variables are as follows: 

H horizontal distance from mesh point to blade, which is hg or h 4 of fig. 17 

K real code variable changed to 1 by INTPL1 

The internal variables are as follows: 

A XI - XO (see fig. 16) 

C X2 - XO (see fig. 16) 

H distance from ZO to approximate root ZINT determined by linear interpolation 

(see fig. 16) 
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HAN length of interval in which root has been determined to lie (see fig. 16) 

IAPN index of vertical mesh line to right of interval 

IAPNM1 index of vertical mesh line to left of interval 

SRW integer variable in COMMON used to obtain output useful in debugging; when 

SRW = 19, values of IA, N, HA, XI, X2, XO, and ZO are printed to start, 
followed by values of H, XI, XI, ZO, and ZBASE for each iteration, and 
final value of H after convergence 

XI 0 -coor dina te of blade computed by BCDCR1 for m = ZINT (see fig. 16) 

XO 0 -coordinate of blade at left end of interval (see fig. 16) 

X2 0-coordinate of blade at right end of interval (see fig. 16) 

ZBASE m-coordinate of left end of interval for first iteration 

ZINT m-coordinate determined by linear interpolation (see fig. 16) 

ZO m-coordinate of left end of interval (see fig. 16) 


SUBRU UT IN E INTPL1 ( X, I A , XI ,H , N ,K , Z , X SP , A LI , ALO, NS P ,S LOPE , EM, RI , RO, 
1 CHORD, STGR, PITCH, SURF, HA ,HB,MXBI .MXBO.RMI ,RMG) 

UIMEN S ION XI 10G) * Z ( NSP ) » X SP ( NSP ) tSLQPE(NSP) ,EM { NSP ) 

COMMON SRW 
INTEGER SRW 
REAL K 
K = 1. 

IAPNM1 = IA+N-1 
IAPN = IA+N 
XO = X ( IAPNM1 ) 

X 2 = X ! IAPN I 
HAN = FA 

ZO = FLOAT! IAPNM1-MXBI )*C HURD /F LOA T ! MXBO-MXB I ) 

IF! IA.ED.MXBO) ZG = CHORD-HA 
H = HA 

IF! IAPN.EQ.MXBO .AND . SURF . LE . 0. . AND. N. NE.O) X2=X2-R0 
IF! SRW. EU. 19) WR I TE ( 6 , 10 10 ) I A , N ,H A , XI , X2 , XO , ZO 
IF! ZO . LT.C.. OR. ZO.GT. {CHORD-. 001*HA) ) RETURN 
ZBASE = ZC 

IF! ABS! Xl-XC ) .GT. ! .001*H8) J GO TO 10 
IF! IA . EO .MXB 1+ 1 ) GO TO 15 
A = 0. 

C = H 
GO TO 20 
15 H = 2 . *R I 
A = -1. 

GO TO 25 
10 A = X1-X0 
C = X2-X0 
2C h = A / C*HAN 
25 CONTINUE 

IF! SRW. EQ .19) WR I TE ! 6 , 1 G20 ) H , XI ,X1 ,Z0 , ZBASE 
ZINT = ZO+H 

CALL BLDCR1 IZ ,XSP, SLOPE ,EM,NSP,RI »ALI , R C,AL0, CHORD, STGR, P ITCH, 
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1 SURF.X I, 1,Z INT.MXB I ,DXDZ,RMI ,RMO) 

I Ft Ati StXI— Xl).LE.tHB*.001) ) GO TO 40 

IFt A*t XI-XI) .LT.O. ) GO TO 30 

HAN = F 

X2 = X I 

GOTO 10 

30 HAN = FAN-H 
XO = X I 
ZO = ZC+H 
GO TO 10 

40 H = Z0+H-Z8ASE 

IFtN.EQ.O ) H = HA-H 

IFt SRW .EQ.19) WR I TE t 6 » 1 GCO ) H 

RETURN 

1000 FORMAT tlX,22HH AS COMPUTED BY INTPL / 1 5X ,i>E 16 .8 ) ) 

1010 FURMAT t IX » 4HIA =,I4,5X,3HN = ,I4 t 5X,4HHA = ,E14.6 ,5 X.4HX1 =,E14.6, 

1 4X.4HX2 =,E14.6,4X,4HX0 = ,£ 14. 6 ,4 X ,4HZ0 =,E14.6) 

1020 FORMAT tlX,3hH =, fc 14. 6 , 5 X , 4H XI * , E 1 4. 6 ,5 X ,4HX1 = ,E14.6 ,5X ,4HZ0 =, 
1 E14.6, 5X.7HZBASE =,E14.6i 

END 


Subroutine VELOC 


The partial derivatives du/dm along each horizontal mesh line are calculated by 
VELOC. This subroutine is described in reference 5. There are no changes in the sub- 
routine, however WX refers here to 0u/9m rather than W^. 


SUBROUTINE VELOC t U ,ZPL , IA1 ,1 AN ,H A1 ,HA K ,1 , 1 B * MX , NXN ,NBB0 , JU , JL, 

1 KK 1 ,KN f USP i RX f RXU.ZXU.RXL t ZXL f NU ,NL ,AAA) 

CIMENS ION UtNXN ) ,ZPLt ICC) , KX (NXN) ,RXL (100) , ZXUtlOO ) ,ZXLt 100 ), 

1 ZAt 100), USPt 100) ,WXSP tlOO) ,WXL(1Q0> ,NU(100) ,NLtlOO) , AAA t 100) 
C KK 1 OR KN = It LOR ER SURFACE 
C KK1 OR KN = 0, UPPER SURFACE 
II = I 
IA2 = IAl + l 
IANM1 = IAN-1 
Z At I A 1 ) = ZPLt IA2J-HA1 
ZAt IAN ) = ZPLt IANM1)+HAN 
NSP = IAN- IA 1 + 1 
00 1C 1A = IA2 » IANM 1 
USPt IA ) = UU1) 

1 1= 1 1+NL t IA )-NUt IA + D+l 
1C ZAt IA I = Z PL t IA ) 

II = NXN+IB-NBBO 
USPt IA 1) = C .0 
USPt IAN) = 0 . 

IFtIAl.EQ.1) USPtl) = UtIB+1) 

IFtKKl.NE.C) USPtlAl) =1.0 
IFt IAN .EQ .MX ) USPt IAN ) = Util) 

IFtKN.NE.C) LSPtlAN) = 1. 

CALL SPLINE t ZA t IA 1 ) , USP 1 1 A 1 ) ,N SP , W X S F 1 1 A1 ) , A A A) 

II = I 

DU 20 IA= IA2 t IANM1 
WXt II ) = WXSPt IA ) 
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2 C II = 1 1+NLI IA )-NU( IA + 1)+1 
C TAKE CARE OF FIRST POINT 
IF( IA 1 .NE .1 ) GO TO 30 
WXI 18 + I ) = WXSPI IA1) 

GO TO fC 

3 C IFIKK1.NE.C) GO TO 40 
JU = Jl + 1 

WXUIJU ) = WXSPIIAl) 

ZXUI JU J = ZA ( IA 1 ) 

GO TO tO 
4C JL = JL+1 

WXLIJL ) = -WXSP( IA1) 

ZXLIJL ) = ZA( I A 1 1 
C TAKE CARE OF LAST POINT 
50 IF( IAN .NE.MX ) GO TO 60 
II = NXN+IB-NBBO 
WXIIl) = WXSPI IAN ) 

RETURN 

6C IFIKN.NE.CI GO TO 7C 
JU = JC + 1 

WXUIJU] = WXSPIIAN) 

ZXUIJU ) = ZAI IAN) 

RETURN 

7 C JL = JL + 1 

WXLIJL } = WXSP I IAN) 

ZXLIJL ) = ZA I IAN ) 

RETURN 

ENO 


Subroutines SPLINE, SPLN22, SPLINT, SORTXY 


These subroutines are all described in reference 5. 


SUBROUTINE SPLINE I X , Y , N , SLOPE ,E MJ 
DIMENSION XI N), YIN), EM IN) , SLOPE IN) 

COMMON Q/ BOX /SI IOC) f A 1100) ,B (100) ,C (100) ,F (100) ,W(100 ) ,SBI 100), 
1GI 100) 

INTEGER Q 
00 10 1 = 2, N 
10 SI I ) = X I I)-XI 1-1) 

NO=N- 1 

IFINO.LT. 2) GO TO 25 
DO 20 1=2, NO 
A I i)=S(I)/i. 

B( I ) = I SI I ) + S( 1 + 1) )/3. 

c 1 1 ) =s i i+n/ 6 . 

20 FI I ) = ( Y( I + l)-Y( I ) )/S(I + l)-( YU l-YII-1) )/S(I) 

25 AIN ) =-I.C 
BI 1 ) = 1 - 
BIN )=1 . 

C( i)=- 1.0 
F( 1) = 0 . 

F I N ) = 0 . 
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W( 1 ) = B l 1) 

SB( l)=Ct 1 )/W( 1) 

G( 1) = C . 

00 30 1 = 2, N 

wi u=b( i j — A t ij*sb ( i— i » 

SB( I)=C( I )/W( I ) 

30 G( I ) — t Ft I )— A (I)*G(I— 1))/W(I) 

EMI N ) = G(N ) 

DO 40 1=2, N 
K=N + 1- I 

40 EM(K)=C(K )-SB(K)*£MIK + l> 

SLUPE( 1J=-S( 2>/6.*(2.*EM(l)+£M(2) ) + (Y(2)-Y(l) )/S(2> 

00 50 I =2, N 

50 SLOPE! l) = S(I)/6.*(2.*EHII)+EH(I-im-(Y(I)-YII-l))/SU) 

IF (U.EQ.13) WRITE (6,100) N ,( X( I ) , Y (I ) ,SLGPE { I) ,EM { I ) , 1= l , N) 

100 FORMAT (2X15HNO. OF POINTS =I3/10X5HX 13X5HY 15X5HSLGPE15X5H 
1EM / I 4F 20 • 8 i ) 

RETURN 

END 


SUBROUTINE SPLN22 ( X , Y , Y IP , YNP , l\ , StCPE ,EM) 

OIMENS BN X( 50) ,Y( 50), SI 50) ,A(5C) ,B (50) ,C(50) ,F(50) , W ( 30 ) , S B ( 50 ) , 
1 G ( 50), £M( £ C ) , SLOPE ( 50 ) 

COMMON Q 
INTEGER Q 
00 10 1 = 2, N 

1C S( I )=X( I ) — X ( I- 1 ) 

NO=N- 1 

IF(NJ.LT.2) GO TO 25 
CO 20 1=2, NO 
A( I) = S( I ) / 6 . 

8( I )=( S( I )+S( I+l))/3. 

C( I ) = S( 1+1 )/6. 

20 F( I ) = ( Y( I + 1J-YI I ) )/S(I + l>-( Y(I l-Y(I-l) )/S(i) 

25 AIN) = S( N )/£ • 

B( l)=S(2)/2. 

BIN) = S(N )/2. 

C( l)=S(2)/6. 

F( 1)=( Y(2)-Y( 1))/S(2)— YIP 
FIN) = YNP-I YIN )-Y(N- 1) ) /SIN) 

W( 1)=B( 1) 

SB! 1 » = C l 1 )/M( 1 ) 

G( 1)=K 1)/W( 1) 

CU 30 1=2, N 

W I I ) = B( I )-A( I )* SB I I— 1 ) 

SB! I)=CI I J/WII ) 

3C G( 1 1 = 1 FI I )-A III*G(I— 1))/W(I) 

EM(N) = C(N ) 

CO 40 1 = 2, N 
K =N + 1— 1 
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4C £M(K)=G(K )— SB ( K }*EM ( K+ 1) 

SLUP£( 1)=-S( 2)/6.*(2.*EM(l )+EM(2) )+ ( Y (2)-Y(l) )/S(2) 

CO 50 1=2, N 

50 SLOPE! 11= S( I )/6.*(2.*EM(I )+£MI I— 1) )+(Y(I)— Y(I— 1) )/S ( I ) 

IF (U.EQ.J8) WRITE (6,100) N,(>!I) ,Y(i) , SLOPE (I) , EMI I ) , 1= 1 , N ) 
RETURN 

100 FORMAT (2X15FNO. OF POINTS =I3/ICX5HZ 15X5HX 10X10FDERIVAT IV 
1E10X10F2ND UERIV./I4G2C. 8) ) 

END 


SUBROUTINE SPLINT ( X , Y ,N ,Z ,MAX , YI NT ,D YD X) 

DIMENSION X(N),Y(N) ,Z(MAX) , YI NT ( MA X) ,D YDX { MAX) 

COMMON Q/80X/S! 50) ,A (50 ,B ( 50) ,C (50) ,F (50) , W (50 ) ,S B (50 ) , G( 50 ) , 
1£M( 400 ) 

INTEGER Q 
III = G 
DO 10 1=2, N 
10 S( I )=X( I ) -X( 1-1) 

NO=N— 1 

IFIN0.LT.2) GO TO 25 
DO 20 1=2, NO 

AI I ) = S ( I )/6.C 
B( I ) = ( S( I ) +S( I* 1 ) )/3. C 
C( I) = S( 1*11/6.0 

20 F( 1) = ( Y( I + I )- Y ( I ) )/S( 1+1)— (Y(I )— Y(I-l) )/S(I) 

25 A(N) =-.5 
B( 1)=1 .0 
B( N ) = 1 .0 
C< 1)=- .5 
F( 1 ) = 0 .0 
F(N ) = 0.0 
W ( 1 )=B( 1) 

SB{ 1)=C( 1 )/W( 1) 

G( 1) = 0 .0 
CO 30 1=2, N 

W ( I ) = B ( I )~ A ( I)*SB( 1-1) 

SB( I)=C( I )/W( I) 

30 G( I ) = ( F( I )-A( I)*G(I-1) )/W(I ) 

£M( N ) = G(N ) 

DO 40 1=2 ,N 
K =N + 1- I 

40 EMIK )=G(K )-SB(K)*EM(K+l) 

DO 90 1=1, MAX 
K = 2 

I F( Z ( I l-X ( 1) ) 6C, 50, 7 C 
50 Y IN T ( I )=Y ( 1 ) 

GO TO 67 

6C IF(Z( I ).GE.l 1.1*X(1)-.1*X(2))) GO TC 85 
WHITE (6,1000) Z(I) 

Q = 16 
GO TO £5 
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1000 FORMAT I17H OUT OF RANGE Z =F1C. 6) 

6 5 K=N 

IFIZI I ).LE.l 1.1*XIN)-.1*XIN-1) ) ) GO TC 85 
WRITE (6, 1000) Z { I ) 

Q = 16 
GO TO 65 

70 IF(Z( I)-X(KI) 65, 75,80 
75 Y INTI I )=Y ( K ) 

GO TO 6 7 
80 K=K+1 

IF(K-N) 7C,7C,65 

85 Y INTI I ) = EMIK- 1)*IXIK)-Zl I ) )**3/6. /SIK) + EMIK) *IZI I l-XIK-l ))**3/6. 
l/SIK)+IYIK)/S(K)-EMlK)*SlK)/6. )* IZI IJ-XIK-1 ) )♦ l Y IK-1 )/ S l K) -EMI K-l ) 
2+SIK )/£. )*l XIK )-Z( I ) ) 

8 7 GY OX I I )=- EMI K— 1 )* I XIK)-Z(I ) )**2/2. G/S I K) + EM I K) * IX l K-l ) -Z l I ) )**2/2. 

10/SIK ) +1 YIK )— YIK-1) J/SIK)— IEMIK)— EMI K— 1 ) ) *SIK)/6.0 
90 CONTINUE 

MX A = MAXC IN ,MAX ) 

I FI U.EC.16) WR ITE (6, 1010) N ,MA X , I X II ) , Y II ) , Z 1 1 ) , Y I NT I I ) ,DY OX I I ) , 

1 1=1, MXA) 

Q = I I I 

1010 FORMAT I2X21HN0. OF POINTS GIVEN =,I3,30H, NO. OF INTERPOLATED POI 
1NTS =, 13, /10X5HX 15X5HY 12 XI 1 H X- 1 NTE RPCL.9 XI 1 HY- 1 NT ERPOL . 

2 8X14HDYDX— INTERPOL. /I5E20. 8) ) 

100 RETURN 

END 


SUBROUTINE SORTXYI X, Y , NPTS ) 
DIMENSION XI 100) » Y( 100) 

100 N=NPTS 
102 NN=N— l 
10 A DO 140 KT=1»NN 
XM IN=X( XT ) 

JAD=KT 

JKL=KT+1 

112 DO 120 JK=JKL,N 

114 IF (XMIN-X(JK)) 120,120,116 

116 XMIN=XIJK) 

118 JAD=JK 
120 CONTINUE 
122 YM IN=YI J AO ) 

XI JAD 1 = XIKT) 

YIJAD)= YIKTI 
XI KT ) = XM IN 
Yl KT ) =YM IN 
140 CONTINUE 
RETURN 
END 


f 
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Lewis Library Subroutines TIME1, ARERR, and DEDERR 


These three subroutines are part of the Lewis Systems Library. TIME1 gives the 
time in clock pulses of 1/60 of a second. To get elapsed time in minutes, the clock must 
be read twice and the difference divided by 3600. TIME1 may be replaced by a user's 
clock reading subroutine, or it may be removed from the program. 

DEDERR and ARERR are error subroutines which return control to the monitor after 
printing out an error message and a trace back giving the location at which the error oc- 
curred. ARERR has the additional feature that it can be overriden by using a "continue” 
control card at the beginning of the deck. The error message is still printed. These sub- 
routines should be replaced by a similar type of error return subroutine at the user’s in- 
stallation. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, November 17, 1967, 

720-03-01-35-22. 



APPENDIX A 


FINITE-DIFFERENCE APPROXIMATION 

An approximate numerical solution for the stream function u can be obtained by 
finite -difference methods. These methods involve first establishing a rectangular grid of 
mesh points in the region as shown in figure 12 (p. 16). Then at each point where the 
value of the stream function is unknown, a finite-difference approximation to equation (1) 
can be written. Adjacent to the boundary, the boundary conditions are included. If there 
are n unknown values, n nonlinear equations are obtained in n unknowns. The equa- 
tions are nonlinear since the coefficients involve the density, which depends on the solu- 
tion. The equations may be solved by an iterative procedure. 

First, the inlet absolute total density is used for determining the coefficients. This 
results in n linear equations. These linear equations may be solved iteratively by 
successive overrelaxation as described in reference 5. There are two major levels of 
iteration in the solution. The inner iteration consists of the iterative solution of n linear 
equations by successive overrelaxation. This solution is an approximate solution of equa- 
tion (1) for the stream function. This approximate solution may be differentiated numeri- 
cally and approximate velocities obtained from equations (2) and (3). The approximate 
velocities are then used to obtain a better approximation to the density at each point, and 
the coefficients of equation (1) are recalculated using new densities. Thus, the solution 
to the nonlinear equation (1) is approached by a sequence of solutions to linear equations. 

2 

h 2 

x h 3 h 4 

3 0 ' 4 

h l 

( 68 ), - 111 
1 r n 


Figure 17. - Notation for adjacent mesh points 
and mesh spaces. 


A typical mesh point with the numbering used to indicate neighboring mesh points is 
shown in figure 17. The value of the stream function or the other variables at 0 is de- 
noted by using the subscript 0, and similarly for the neighboring points. It can be shown 
(ref. 13) that equation (1) can be approximated by 
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where = Tq(A0)j and hg = ^(ASjg (since Tq = r^ = rg). For settingup our equations 
for solution, the coefficients of the u. in equation (Al) must be calculated. This was 
done by expressing equation (Al) as 


where 


u o = Z a i u i + k o 

i=l 


h l h 2 


d 34 


h 3 h 4 


a 0 " a 12 + a 34 


b 12 _ 


p 2 ~ P 1 
P 0 ( h l + h 2> 


^4^4 ” ^3^3 s ^ n 


34 


boP 0 ( h 3 + h 4 ) 


1 a o< h l + h 2 > \ h l 


— + b 


12 


d 12 


i 2 h 

*3 -r,~lr + b 34 


a 0 (h 3 + h 4> \ h 


_ 34 „ 

* 4 “ “ " a S 
a 0 


k o = -^Vo sin "o 


> (A 2) 
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This equation can be used at all interior mesh points, and for mesh points adjacent to the 
blade surfaces BC or FG. 

Along the boundary where the value of u is unknown, the equation will vary. For 
example, along the upstream boundary, 3u/3 77 is known, and a finite-difference approxi- 
mation to (3u/377) in in equation (5) gives 


'tan (3. 


u 0 = u 4 + -4+* 4 r- n 

\aVin 


sr in 


Similarly, along the downstream boundary, 


u 0 = u 3 + h 3 


^ =u 3 -h 3 f?^ 

?n)out \ sr out > 



(A3) 


(A 4) 


U 1 " u l, s - 1 



Figure 18. -Mesh point on line AB. 


For the points along AB and CD equations can be derived by using the periodic 
boundary condition. If the point 0 (fig. 18) is on the boundary between A and B, the 
point 1 is outside the boundary. However, it is known that u^ = u^ g - 1 where the 
point l,s is a distance s above point 1 in the 9 direction, as shown in figure 18. Sub- 
stituting this condition in equation (A2) gives 

4 

u 0 =a l u l,s + Z a i U i- a l +k 0 (A5) 

i=2 
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where the are the same as defined in equation (A2). Of course, equation (A 5) holds 
along CD also. 

The points along GH need not be considered, since they are just 1 greater than the 
corresponding point along AB. The equation for the first mesh line below HG must be 
modified. In this case u 2 = u 2 , + 1, where the point 2, -s is a distance s below 

point 2 in the negative 9 direction, as indicated in figure 19. Substituting this condition in 
equation (A2) gives 


u 0 = a lUj + a 2 u 2) _ g + a 3 u 3 + a 4 u 4 + a 2 + k Q 
Equation (A6) also applies to the first mesh line below FE. 


(A 6) 


H 2 




u 2 “ u 2 -s + 1 



Figure 19. - Mesh point on first line below HG. 


One of equations (A2) to (A6) can be applied to each mesh point for which the stream 
function is unknown in the region of interest, giving the same number of equations as 
there are unknowns. These points where the stream function is unknown will be referred 
to simply as unknown mesh points. 

This system of n equations is represented in matrix form as 

Au = k (A 7) 

T 

where u = (uj, . . . , u n ) is a vector whose components are the unknown values of the 
stream function, A is the coefficient matrix of equations (A2) to (A6), and k = (k 1 , . . . , 
k R ) 1 is the vector whose components are the known constants of equations (A2) to (A6). 

If the mesh size is sufficiently small, the coefficients, a 1 to a 4 in equation (A2) will be 
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all positive (for any given continuous function b and p). In this case, the coefficient 
matrix A is irreducibly diagonally dominant, and there is a unique solution to equa- 
tion (A2) (ref. 13). 

The solution to equation (A2) is obtained by using two levels of iteration. The inner 
iteration consists of solving (A2) using fixed values of p based on the previous inner 
iteration. The inner iteration is successive overrelaxation using an optimum overrelaxa- 
tion factor fi, as described in reference 5 (p. 77). The iterative procedure is given by 


u m+1 

l 



+ 


i-1 


I 


a..uP +1 
i] 1 


i vr +k i 

j=i+l 



for i = 1, 2, . . . , n (A8) 

where is the overrelaxation factor. The a., are the elements of the matrix A, and 
the are the components of the vector k of equation (A7) . The u. are the initial es- 
timates of the m and are obtained from the previous inner iteration. The optimum 
value of can be determined as described in reference 5, appendix B. The optimum 
value of S7 will vary slightly each time the coefficients are corrected; however, the 
change is usually small, and it has been adequate to use the same overrelaxation factor 
for the entire calculation. 
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APPENDIX B 

NUMERICAL TECHNIQUES USED IN PROGRAM 
Calculation of Velocity and Density 

When the stream function u has been calculated, it is possible to then calculate the 

derivatives 3u/3m and 3u/30 by numerical techniques. Then, with equations (2) and (3), 
2 2 2 

and since W = W m + W^, values for pW can be calculated. It is assumed that the 
values of o, X, r, y, c , T? n , and p? n are all fixed and known. Then p, and hence 
pW, is a function of W. The product pW has its maximum value when W-W cr _ If pW 
is less than this maximum value, there are two values of W which will give this value of 
pW, one being subsonic and the other supersonic. It is desired to find the subsonic value 
of W corresponding to the given value of pW. The method used in Newton's method, 
which converges quadratically. 

It is necessary to express pW as a function of W. Since 

V 2 = W 2 + 2oA - (or) 2 


we have 



W 2 + 2oA - (or) 2 

2 c T! 
p in 


With the assumption of isentropic flow 


1 


P - 


p'. 



(Bl) 


(B2) 


hence, 


pW 


Pi 


in 


1 - 


W z + 2pA - (or)' 

2 c t: 
p in 


1 y -1 


w 


(B3) 
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For Newton’s method, the derivative with respect to W is needed, 


2-y 1 


d(pW) w2 fin 

j W 2 + 2a>A - (wr) 2 

y- 1 

+ p! 

'in 

x _ W 2 + 2wX - (wr) 2 

dW yRTjn 

2 c T! 

L pm J 

2 c T! 

L p in J 


Suppose that (pW)^ v is a given value of pW. A first estimate of W is 


(PWW 

\xr - & 1V 

° —p'- 


Then, using Newton's method, 


W n+! = W n + 


(PWlgiv- 

■ P(W n )W n 

d(pW) 


dW 

w=w„ 


n = 0, 1, 2, . . . 


n 


(B4) 


(B5) 


(B6) 


Since the convergence is quadratic only a few iterations are needed and the relative 

change in IV is an excellent measure of the relative error in W . If an estimate for 
n n 

W is available from a previous iteration, then this value is used for Wq instead of 
using equation (B5). The algorithm given by equation (B6) is done by subroutine DENSTY. 


Calculation of X 


The input information for the program determines the value of 
value of (pW)^ n can be calculated by 



The 


<pW)i„ = 


w 

r in sb in cos /3 in 


(B7) 


since there is assumed to be uniform flow across AH. The value of W can be estimated 
by dividing this value of (pW) in by p! . Then X can be estimated by 

X = r in( W in sin ^in + wr in> ( B8 > 
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and from this a better value of p is calculated by 


1 

W 2 + 2cuA - (cur) 2 


Use of this value of p gives a better estimate of the value of W, and then iteration can 
be used with equations (B8) and (B9) until there is a negligible change in p. This cal- 
culation also gives the value of W along AH. These calculations are performed in 
INPUT, after reading all input cards. 


y-i 


(B9) 


Calculation of W cr 

For reference the critical relative velocity W cr is calculated at blade leading and 
trailing edges. This is given by 



_ 2yR T „ 
y+1 


(BIO) 


where 


T f 1 = T 


in 


2cuA - (cur)' 
2 c 


This calculation is performed by COEF after reading the input cards. 


(Bll) 


Calculation of Maximum Value of Mass Flow Parameter pW 


The mass flow parameter pW attains its maximum value when W = W . For ref- 
erence, the maximum values of pW along AH and along DE are computed by the pro- 
gram. The maximum value of pW is calculated by 

1 


(pW) = p? 
^ 'max 'in 


1 - 


W 2 r + 2cuA - (cur) 2 


2 c T! 
p in 


y-i 


w 


cr 


(B12) 
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where W cr is calculated by equations (BIO) and (Bll). 


Calculation of |3 at Leading and Trailing Edges 

If the radius r or streamsheet thickness b is not constant, /3 may change for free- 
stream conditions. At the inlet, the hypothetical freestream angle /3 may be calculated 
by 


tan /3 = ( A - <*> r2 ) p bs 
w 

Equation (B13) may be derived from the following relations for uniform flow in the 
9 direction: 


(B13) 


tan /3 


W, 


w 


m 




- wr 


A = rV 0 


W 


m 


w 

pbrs 


>- 


J 


(B14) 


At the exit rV 0 is constant (/\) and the other relations in equation (B14) hold. Using 
this gives 


tan /3 


tan 0* / p \ w(r* - r 2 )ps 

b* \P* J w 


(B15) 


for the freestream angle (3 where * denotes values at some reference coordinate of 
m = m* . 

Equation (B15) may be used at either inlet or outlet to calculate /3. or /3 t> In the 
program, equation (B13) is used to calculate the freestream angle /3 at the leading edge 
BG, and equation (B15) is used to calculate the freestream angle /3 at the trailing edge 
CF. In the program p in is used for p in equation (B13) and p Qut is used in equa- 
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tion (B15) since there is little change in density in the freestream region. 


Equation for Leading- and Trail i ng-Edge Radii 

The equation for the leading- and trailing-edge radii is needed. If the radius r were 
constant, 


(m - m*) 2 + r 2 (d - 0*) 2 = R 2 (B16) 

where R is the leading- or trailing-edge radius and (m* , 0*) are the coordinates of the 
center of the radius. Since r changes by a relatively small amount on this circle, it 
was deemed adequate to use this equation with r taken at the leading or trailing edge. 
Equation (B16) is used by the program to calculate coordinates on the leading- and 
trailing-edge radii. It is also used to calculate the points of tangency to the spline curves 
describing the rest of the blade surfaces, and to calculate slopes on the leading- and 
trailing-edge radii. 


Calculation of Surface Length 


It is often desired to plot the velocities as a function of blade surface length. For 

convenience, the approximate blade surface length is calculated by the program. The 

calculation is based on straight line distances between each vertical grid line on the blade 

th 

surface. If h A is the spacing between vertical grid lines, r. the radius at the i ver- 
tical grid line, and 6. the coordinate of the i vertical grid line, the surface length S 
to the n grid line is approximately 



(BIT) 


This may be in error near the leading or trailing edge, but is quite accurate over most of 
the blade surfaces. 
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